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Summary 

Developing  state  bound  estimation  algorithms  for  nonlinear  systems  has  been 
of  high  importance  in  robustness  analysis  of  dynamic  systems.  For  many  cases, 
Monte-Carlo  simulation  might  be  the  only  tool  to  estimate  these  bounds  for  a 
general  type  of  nonlinear  systems.  The  required  number  of  simulations  for  a 
tight  bound,  however,  would  be  very  large  and  it  might  be  impossible  to  com¬ 
plete  within  a  given  computation  time,  /i-formulation  for  state  bounds  trans¬ 
forms  the  bound  estimation  problem  to  a  singularity  problem  and  the  singular 
problem  is  solved  using  a  randomised  optimisation  approach.  The  performance 
of  the  algorithms  are  demonstrated  by  multi-dimensional  Rosenbrock  function; 
simple  discrete  system;  large-scale  biological  system;  hybrid  system;  and  nav¬ 
igation  error  propagation  for  underwater  vehicle.  For  a  given  error  tolerance 
of  the  bounds,  a  formula  to  calculate  the  required  number  of  sampling  in  the 
algorithms  is  provided.  Because  of  the  inherent  complexity  of  general  nonlinear 
optimisation  problems,  the  required  sampling  number  increases  very  fast  as  the 
problem  dimension  increases.  The  suggested  algorithms  would  produce,  how¬ 
ever,  tighter  estimation  faster  than  random  blind  search.  In  addition,  exploiting 
parallel  computation  architecture  the  suggested  algorithms  could  be  the  solu¬ 
tion  of  real-time  robustness  analysis  in  future. 

The  following  two  conference  papers  are  accepted  to  19th  IFAC  World  Congress 
2014,  24th-  29th  August  2014,  Cape  Town,  South  Africa: 

1.  Kim,  J.,  Kishida,  M.,  and  Bates,  D.G.,  ’’State  bounds  estimation  for  non¬ 
linear  systems  using  mu-analysis” 

2.  Kim,  J.,  Park,  Y.,  Lee,  S.,  and  Lee,  Y.K.,  ’’Underwater  glider  navigation 
error  compensation  using  sea  current  data” 
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Chapter  1 


Introduction 


Robustness  is  the  capability  of  systems  maintaining  their  functions  while  in¬ 
ternal  and  external  disturbances  produce  adverse  effects  on  the  performance. 
Robustness  analysis  is  an  indispensable  step  in  designing  engineering  systems 
[1,  2,  3].  It  is  also  considered  as  one  of  the  most  important  aspects  in  analysing 
biological  systems  [4,  5,  6,  7,  8,  9,  10].  Systematic  approach  to  model  many 
interactions  and  analyse  noisy  measurement  data  is  one  of  the  highly  preferable 
ways  in  improving  understanding  of  complex  systems  [11,  12]. 

Recently,  some  of  biological  system  models  are  described  with  hundreds  of 
states  and  parameters  [13,  14].  And,  such  large-scale  systems  provide  infor¬ 
mation  that  was  not  available  with  small  scale  models.  Similarly,  engineering 
systems  have  been  evolving  to  be  more  and  more  complex  in  order  to  accomplish 
improved  performance  while  attaining  stronger  robustness  towards  various  un¬ 
certain  environment.  This  tendency  will  continue  in  future  for  satisfying  much 
demanding  design  specifications  and  it  is  important,  hence,  to  have  efficient 
numerical  methods  to  analyse  large-scale  systems. 

Structured  singular  value  or  /i-analysis  has  been  one  of  the  most  successful 
tools  for  robustness  evaluation  [15,  16,  17,  18].  Even  though  the  computational 
complexity  of  /x-analysis  was  acknowledged  at  the  beginning  state  of  the  concept 
introduced,  it  has  been  successfully  used  for  many  practical  systems  design  and 
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analysis. 

The  computational  complexity  of  /x-analysis  is  proved  in  [19]  by  transform¬ 
ing  /x  problem  to  a  corresponding  optimisation  problem  that  is  known  to  be 
Non-deterministic  Polynomial-time  hard  (NP-hard).  As  far  as  NPy^P,  the  com¬ 
putational  complexity  is  the  fundamental  obstacle  that  cannot  be  overcome  by 
any  algorithms  without  introducing  conservatism  in  design  and  analysis.  Practi¬ 
cally,  however,  /x-analysis  algorithm  produced  many  useful  results.  In  addition, 
recently,  it  is  further  extended  to  solve  some  class  of  optimisation  problems  [20] . 
This  is  an  inverse  interpretation  of  the  formulation  in  [19]. 

The  application  in  [20]  is  calculating  state  bounds  for  polynomial  nonlinear 
discrete  systems,  where  initial  states  and  parameters  in  the  systems  are  de¬ 
scribed  by  uncertain  bounds.  Then,  the  state  maximum  and  minimum  bounds 
are  calculated  using  /i  upper  and  lower  bounds  algorithms.  As  long  as  the  non¬ 
linearity  appears  in  polynomial,  the  uncertainties  can  be  decoupled  from  the 
known  parts  and  the  system  can  be  described  in  Linear  Fractional  Transforma¬ 
tion  (LFT)  [19,  21].  Once  LFT  form  is  obtained,  some  powerful  numerical  tools 
that  provide  the  upper  and  the  lower  bounds  of  /.t  can  be  used,  for  example,  /z 
toolbox  or  Robust  Control  toolbox  in  MATLAB  [1]. 

Conservatism  of  the  calculated  bounds  and  the  requirement  for  the  uncer¬ 
tainty  structures  for  enabling  LFT  format  are  two  main  obstacles  for  the  state 
bounds  algorithm  to  be  further  extended  its  applications,  yx-analysis  problem 
is  interpreted  in  geometrical  point  of  view  in  [22]  and  it  enables  one  to  use 
random  sampling  approach  to  obtain  the  bounds.  Later,  it  is  further  lifted  the 
requirement  of  LFT-transformation,  i.e.  LFT-free  yx-analysis  [23].  It  is,  hence,  a 
natural  fusion  that  combining  state  bounds  estimation  and  LFT-free  /z-analysis. 

In  this  research  calculating  the  bounds  of  a  polynomial  function  using  the 
skewed  /x-analysis  framework  presented  in  [20]  is  extended  to  a  general  type  of 
nonlinear  systems  including  discontinuous  and  non-smooth  nonlinear  functions 
using  a  random  sampling  method  [23]. 

This  report  is  organised  as  follows:  Firstly,  state  bounds  estimation  prob- 
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lem  is  formulated  as  LFT-free  ^-analysis  and  state  worst-bounds  algorithms  are 
presented;  Secondly,  the  performance  of  the  algorithms  is  demonstrated  using 
various  examples:  hello  world  example,  oscillatory  discrete  system,  large-scale 
biological  system,  hybrid  system,  and  underwater  glider  navigation  error  prop¬ 
agation,  where  some  of  the  examples  are  used  the  algorithm  parallelised  to  run 
on  Graphical  Processing  Unit  (GPU);  Finally,  the  conclusions  are  presented. 
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Chapter  2 

Methods,  Assumptions,  and 
Procedures 

2.1  Bounds  estimation  in  LFT-formulation 

A  nonlinear  dynamical  system  is  given  by 

x  =  f(x,p),  (2.1) 

where  Xo  =  a;(0)  is  the  initial  condition  at  time  t  =  0,  x  is  the  derivative  of  x 
by  time,  t  ,  x  is  the  state  vector,  which  is  an  element  of  the  n^-dimensional  real 
number  set,  R"*,  nx  is  a  positive  integer,  p  is  the  uncertain  parameters,  which 
is  an  element  of  rip-dimensional  real  number  set,  R™p,  np  is  a  positive  integer, 
and  f(x,p)  is  a  nonlinear  function  in  x  and  p ,  which  might  have  discontinuous 
and/or  non-smooth  parts.  Many  dynamic  systems  will  fall  into  this  category, 
for  example,  sun  tracking  controller  for  UKube-1  using  magnetic  torquer  [24], 
spacecraft  attitude  control  using  thruster  [25],  auto-tuning  control  mechanism 
for  industrial  applications  [26]. 

For  the  well-posedness  of  the  problem,  f(x,p)  is  assumed  to  have  a  unique 
solution  for  the  nonlinear  differential  equation.  For  a  chosen  positive  real  num- 
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ber,  At,  a  transition  function  ip(, )  is  defined  to  satisfy  the  following: 

xk+1  =if>(xk,p)  (2.2) 

where  xk  =  x(kAi)  for  fe  =  0,1,2,....  The  transition  function  can  be  obtained 
by  a  standard  numerical  integrator  such  as  Runge-Kutta  method. 

The  main  problem  is  finding  the  worst-bounds  for  the  maximum  and  the 
minimum  of  a  function,  < f>(xk ),  for  the  given  bounds  for  Xq  and  p  as  follows: 

Problem  2.1.1  (Bounds  estimation)  For  a  real-valued  scalar  function,  (j>(xk,p), 
obtain  </>  .  ,  4>  and  ^max  in  the  following  inequalities: 

— min  — m  3.x 

$min  <±<  </>  min,  (2.3a) 

—max  —  ^  —  ^maxj  (2.3b) 

where  <f>  and  4>  are  the  unknown  true  bounds  satisfying  f>  <  <fi(xk,p)  <  <j)  for  all 
Xq  <  Xo  <  Xq,  and  p  <  p  <  p.  That  is,  <f>  is  the  minimum  upper  bound  and  f>  is 
the  maximum  lower  bound. 

Remark  2.1.2  The  real-valued  scalar  function,  <j>(xk,p),  could  be  the  position 
bounds  of  a  group  of  debris  in  space,  the  miss  distance  bounds  of  missile,  or  the 
navigation  error  bounds  of  mobile  robot. 

In  the  above,  all  the  upper  and  the  lower  bounds  for  <f>(xk,p),  Xo  and  p  are 
assumed  to  be  finite.  Systems  having  a  finite  escape  time  to  either  positive  or 
negative  infinity  are  excluded  in  this  study 

Let  <p(xk,p)  =  Xk  in  Problem  2.1.1.  Then,  the  problem  becomes  state  bound 
estimation  problem  for  the  give  uncertain  ranges  of  the  initial  condition  and  the 
parameters  in  the  system.  For  the  case  that  f>{xk,p)  is  a  polynomial  function  in 
Xk  and  p,  it  can  be  transformed  into  LFT  form  and  existing  p  bounds  algorithms 
can  be  directly  used  to  obtain  the  bounds.  See  the  below  example. 
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Example  2.1.3  (Special  Case)  Find  the  polynomial  expression  for  the  state 
bounds  at  t  =  2,  i.e.  x(2),  of  the  following  dynamic  system: 


x  =  —(1  +  p)  x 


(2.4) 


where  0  <  x(0)  <  1  and  ——  <  p  < 
sol)  The  transition  function  is  given  by 

</>[x(t),p]  =  x(t)  =  ip[x(t),p\  =  e_<-1+p^a;(0) 


At  t  =  2, 


4>{x{2 ),p\  =  e  (1+p^2x(0)  =  e  2e  2pa;(0) 

Using  the  Taylor  series  expansion,  the  uncertain  exponential  function  is  given 

by 

<j>[x{2),p]  =e-2  ^1-2 P+  \  -  ^(0) 

As  \p\  <  0.5,  ignore  the  higher  order  terms  in  p  as  follows: 

<f>[x(2),p]  «  e-2  (l  -  2p  +p2)  x(0)  (2.5) 

Let  the  nominal  values  for  x(0)  and  p  be  1/2  and  0,  which  are  the  midpoints  of 
the  given  uncertain  boundaries,  respectively.  Using  the  midpoints,  each  can  be 
written  as  follows: 


x(0)  =  \  +  \&* 
P  =  0  +  -Sp 


where  —  1  <  Sx  <  1  and  —1  <  Sp  <  1. 


(2.6a) 

(2.6b) 
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Substituting  (2.6)  into  (2.5) 


ct)[x(2),p\  «  e  2  ^1  -  Sp  +  ^  j  ^  (2.7) 

Hence,  this  is  now  in  a  polynomial  function  and  it  can  be  transformed  into 
LFT-form.  Details  about  the  transformation  of  the  special  case  can  be  found  in 

□ 

On  the  other  hand,  the  Problem  2.1.1  is  not  possible  to  transform  into  a 
polynomial  expression,  in  general,  without  introducing  some  approximation  er¬ 
ror,  which  might  cause  significant  conservatism  in  the  estimated  bounds.  In  the 
following,  the  same  idea  obtaining  the  state  bounds  using  /^-formulation  in  [20] 
is  extended  without  introducing  any  polynomial  approximation. 


Let  Xq  and  p  in  Problem  2.1.1  be  the  following  form: 


x0  —  xc  -I-  (2.8a) 

P  =  Pc  +  WPSP ,  (2.8b) 


where 


Wx 


diag 


(2.9a) 


Wp  =  diag 


wPl 


wP2 


(2.9b) 


diag[. . .]  is  the  diagonal  matrix  whose  diagonal  terms  are  given  in  the  bracket, 
wXi  =  (So  —  x0)/2,  wp.  =  (p  —  p)/ 2,  xc  and  pc  are  defined  such  that 


-  1  <  SXi  <  1  (2.10a) 

-  1  <  SPj  <  1  (2.10b) 


for  i  =  1, 2, . . . ,  nx  and  j  =  1,2,...,  np. 
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Figure  2.1:  Pseudo-LFT  form 


Example  2.1.4  In  the  system  shown  in  Example  2.1.3,  the  original  function 
without  the  approximation  is 


(j>[x(2),p]  =  e  2e  2px0 


where  p  =  pc  +  wpSp,  x0  =  xc  +  wx5x,  pc  =  0,  wp  =  1/2,  xc  =  1/2,  wx  =  1/2, 
and  —1  <  6P  <  1,  —1  <  Sx  <  1. 


Define 


A (fk{Sx,Sp)  :=  cj)k(xo,p)  -  (f)k{xc,pc )  (2.11) 


where 


4>k{-,p)  ■■=  (f[  itoifoifo  ...oil)(-,p),p]  (2.12) 

fc-times 


and 


^(■,P)  =  =  0{V#(-,P),P],P}  (2-13) 

Pseudo-LFT  form  as  shown  in  Figure  2.1  is  to  be  constructed.  It  is  called 
pseudo-LFT  as  it  is  in  the  LFT  format  but  it  can  be  only  evaluated  for  a  fixed 
Sx  and  Sp.  In  the  standard  LFT-formulation,  A <j)k  is  a  constant  matrix  with  a 
structure.  On  the  other  hand,  in  the  pseudo-LFT  formulation,  it  is  a  varying 

11 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


vector  evaluated  on  a  specific  combination  of  6X  and  Sp. 


From  the  pseudo-LFT  shown  in  Figure  2.1  and  the  equivalency  between  p- 
bounds  and  the  optimisation  problem  shown  in  [19],  the  maximum  of  \<j>k(xo,p)\ 
is  bounded  above  as 


max  \(f>k(xo,p)\  <  (2.14) 

tv 

where  k*  is  the  minimum  k  satisfying  the  following  singular  condition: 

\I2-NA\=0,  (2.15) 

|  •  |  is  the  determinant  of  matrix,  /2  is  the  2x2  identify  matrix, 

0  1 

N  = 

l  4>k{xc,pc) 

A  <hk  0 

A  = 

0  kSc 

|5C|  =  |<5r  +  Sij\  <  1,  5r  and  Si  are  the  real  numbers  whose  magnitude  is  less 
than  or  equal  to  1,  and  j  =  \J  ~ 1.  Note  that  k  is  a  positive  number  satisfying 
the  singular  condition. 


(2.16a) 

(2.16b) 


The  singularity  condition  is  expanded 


1 

c 

se 

o 

i _ 

1  —kSc 

II 

<r 

i 

I-? 

A(j)k  K(j)k(xc,pc)5c 

-A cj)k  1  -  K(/)k(xc,pc)Sc 

=  1  -  Kcj)k{xc,pc)Sc  -  nA(j)k5c 

=  {l  -  K  \4>k{xc,pc)  +  A(/)k]  SR}  -  K  \(f>k(xClpc)  +  A 4>k]  Sij  =  0,  (2.17) 

where  the  real  part  and  the  imaginary  part  must  be  equal  to  zero  at  the  same 
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time  for  the  matrix  to  be  singular.  The  real  and  the  imaginary  parts  are  equal 
to  zero  as  follows: 

-  NA\)  =  1-K  [/( xc,Pc )  +  A0fc]  6r  =  0,  (2.18a) 

3(| h  ~  NA\)  =  k  [<j>k{xc,pc)  +  A0fe]  ST  =  0,  (2.18b) 

where  Sft(-)  and  3(-)  are  the  real  part  and  the  imaginary  part  of  the  argument, 
respectively.  Notice  that  (j)k(xc,pc)  +  A f>k  is  equal  to  4>k{xo ,p)  by  the  definition 
and  it  is  not  equal  to  zero,  in  general. 

Hence,  Sj  must  be  equal  to  zero  for  the  singularity  condition.  Therefore,  the 
following  is  the  only  case  that  both  the  real  and  the  imaginary  parts  are  equal 
to  zero:  1  —  k  \cj)k(xc,pc)  +  A^fe]  Sr  =  0  and  <5/  =  0,  i.e.,  the  imaginary  value  of 
Sc  is  always  equal  to  zero. 


Take  the  absolute  value  of  (2.18a)  after  moving  the  second  term  into  the 
right  hand  side  of  the  equation  as  follows:  1  =  |k  <fik(xo ,p)Sr  |,  and  it  becomes 


1 

\(j)k(x0,p)SR  |' 


(2.19) 


Finally,  the  singular  problem  for  the  bound  estimation  is  formulated: 


Problem  2.1.5  (Singular  problem)  Find  the  minimum  k,  i.e. 
the  following: 

K*  =  1  =  1 

max  \<t>k(xo,p)6ft\  max  |^fe(a:o,p)| 


k*  ,  satisfying 

(2.20) 


where  |<5r|  <  1. 


This  looks  trivial  and  it  does  not  seem  to  help  to  find  the  bounds  for  4>k(xo  ,p). 
In  the  next  section,  bounds  algorithms  using  a  random  sampling  approach  will 
be  presented. 


Example  2.1.6  For  the  system  shown  in  Example  2.1.4,  the  singular  problem 
is  given  by  n*  =  1/  max  \e~2e~2pXo\ .  r-i 
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Figure  2.2:  Sign  changes  along  the  uncertain  box  boundary:  as  At  is  smaller  than 
At*,  all  signs  of  g(5x ,  SP:k)  for  the  samples  at  the  boundary  are  positive;  on  the 
other  hand,  when  R  is  greater  than  At*  ,  there  are  red  box  samples  and  blue  circle 
samples  whose  sign  of  g(5x,6p,R)  is  positive  and  negative,  respectively. 

2.2  Worst  Bounds  Estimation  Algorithms 


Define 


g(6x,6p,K)  :=  1  -  At|</>fc (x0,p)|  =  1  -  At \(f>k{xc  +  Wx5x,pc  +  WP5P) |  (2.21) 

where  —r<5x<r,—r<6p<r,K>  0,  and  r  £  (0,1].  Note  that  g(-,-,-) 
could  be  discontinuous  as  <pk{-,-)  could  be  discontinuous.  As  shown  in  Figure 
2.2,  there  are  three  cases: 

•  for  k  <  k*  ,  having  one  type  of  samples,  whose  signs  are  all  positive, 

•  for  At  >  k*  having  two  types  of  samples,  whose  signs  are  positive  or  nega¬ 
tive, 

•  for  At  =  At*  having  at  least  one  singular  point,  where  g(Sx,  5p,  At)  is  equal 
to  zero. 
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Hence,  whenever  for  a  fixed  k,  if  two  sign  combinations  are  found,  it  becomes 
the  upper  bound  for  k*  ,  i.e. ,  R.  Similarly,  if  only  positive  signs  are  found,  it 
becomes  the  lower  bound,  k. 

Example  2.2.1  For  the  system  shown  in  Example  2.1.6,  the  function,  g(Sx,  Sp,k), 
is  given  by 


g{6x,5p,K)  :=  1  k 


□ 


The  bound  estimation  problem  is  given  by 

Problem  2.2.2  (k* -bounds  estimation)  For  the  following  real-valued  scalar  func¬ 
tion: 


g(6x,8p,n)  =  1  -  k  1 4>k{xc  +  WxSx,pc  +  Wp5p)\ ,  (2.22) 

obtain  the  upper  (R)  and  lower  (k)  bounds  as  follows: 

k<k*<R  (2.23) 

Remark  2.2.3  The  upper  bound,  R,  is  a  deterministic  bound  as  we  found  the 
negative  sign  but  k  is  probabilistic  as  it  has  always  some  danger  to  be  failed 
depending  on  the  number  of  samples  checked  on  the  boundary.  A  safety  factor 
could  be  introduced  in  order  to  reduce  the  probability  for  the  upper  bound  to  be 
failed. 

Remark  2.2.4  Once  the  bounds  for  n*  are  obtained,  the  bounds  for  \<pk(-,-)\, 
is  obtained  such  that  <j>  .  (or  <j>  )  is  equal  to  1/R  and  4>min  (or  4>ma,x)  is  equal 

to  1  /  K. 

How  to  resolve  whether  they  are  the  bounds  for  the  minimum  or  the  max¬ 
imum  will  be  presented  later.  In  the  following,  bounds  algorithms  for  a  fixed 
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radius  uncertain  space,  r,  which  is  greater  than  0  and  less  than  or  equal  to  1, 
are  presented  and  an  algorithm  for  sweeping  r  in  (0, 1]  will  be  presented. 

2.2.1  Fixed  uncertain  space  radius:  r  G  (0,1] 

First  part  of  the  algorithm  is  for  calculating  the  pre-K  bound. 

Algorithm  2.2.5  Pre-K  estimation 

1.  Set  N,  the  number  of  samples  along  the  face  of  the  uncertain  box  shown  in 
Figure  2.2,  where  the  size  of  the  box  is  r  measured  in  terms  of  oo  -norm. 

2.  Set  the  initial  boundary  for  k  such  that  e  <  k  <  E,  where  e  could  be  the 
smallest  positive  number  and  E  could  be  the  largest  positive  number  that 
can  be  expressed  in  the  computer. 

3.  Set  the  tolerance,  e,  for  the  magnitude  of  the  interval,  [e,  E\,  i.e.  E  —  e 
Set  initial  guess  of  n,  equal  to  (e  +  E)/2 

5.  For  i  =  1  to  N 

•  Evaluate  g(5x,6p,n)  for  the  given  i-th  sample  of  Sk  and  Sp 

•  If  g{fix,  dp,  k)  <  0,  then  replace  E  by  k  and  break  the  for-loop,  else 
continue  the  for-loop 

6.  If  i  =  N  and  g{8k,  8P,  k)  for  all  samples  are  positive,  then  replace  e  by  k. 

1.  If  E  —  e  is  smaller  than  e,  then  declare  Kp  =  E  and  stop.  Otherwise,  go 
to  step 

Note  that  Kp  from  the  pre-K  estimation  algorithm  does  not  need  to  be  tight 
as  long  as  it  is  smaller  than  k*.  The  reason  to  calculate  Kp  using  the  above 
algorithm  before  calculating  bounds  for  (j)k(x0 ,p)  is  that  the  value  of  f>k(x0,p) 
can  be  positive  and  negative  and  the  definition  of  k*  in  (2.20)  is  given  in  terms 
of  the  absolute  value  of  f>k(xo,p).  On  the  other  hand,  if  the  sign  of  all  possible 
values  of  </>fc  is  either  only  positive  or  negative,  this  step  can  be  skipped. 
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In  order  to  estimate  the  bounds  for  the  maximum  cf>k(xo,p),  the  following  is 
defined: 


g{Sx,dp,K)  :=  1  -  K\(j)k(x0,p)  +  s/kp\  (2.24) 

where  s  is  the  shifting  factor  strictly  greater  than  1.  As  1/kp  >  \&k(xo,p)\ 
can  be  guaranteed  only  in  a  probabilistic  sense,  the  safety  factor  will  make 
sure  the  terms  inside  the  absolute  sign  be  positive.  Then,  the  maximum  of 
\4>k{x0,p)  +  s/np |  occurs  at  f>k(x0,p)  +  s/np. 

Algorithm  2.2.6  (f-bounds  estimation  Algorithm 

1.  Run  the  pre-K  estimation  algorithm  after  replacing  g(5x ,  Sp,  k)  by  g(5x,6p ,  k) 

2.  Set  k  =  e  and  R  =  E 

3.  Declare  <j>  =  1/k—s/kp  and  </>max  =  s//k  —  s/kp,  where  Sf  is  the  safety 

factor  greater  than  1.  The  smaller  Sf  returns  tighter  bounds  with  greater 
danger  to  be  failed.  The  larger  Sf  returns  conservative  bounds  with  less 
probability  to  be  failed. 

Similarly,  the  bounds  for  the  minimum  is  obtained  by  defining 

g(Sx,  Sp,  k)  :=  1  -  n\(j)k(xo,p)  -  s/kp I,  (2.25) 

which  make  sure  the  maximum  occurs  at  the  minimum  of  <pk ,  and  executing  the 
following  algorithm: 

Algorithm  2.2.7  f-bounds  estimation  Algorithm 

1.  Run  the  pre-K  estimation  algorithm  after  replacing  g(5x ,  Sp,  k)  by  g(5x,6p ,  k) 

2.  Set  K  =  e  and  R  =  E 

3.  Declare  4>min  =  1/R  +  s/ Kp  and  ^m;n  =  Sf/n  +  s/ Kp 
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Remark  2.2.8  The  bounds  calculated  by  the  above  algorithms  are  for  a  fixed 
r.  Hence,  the  global  bounds  must  be  obtained  by  finding  the  maximum  bounds 
among  the  ones  for  different  r.  However,  the  lower  bound,  the  worst  scenario 
found  so  far,  is  always  valid  and  the  issue  is  how  close  it  will  be  to  the  true 
bound. 

2.2.2  Sweeping  uncertain  space:  r  e  (0,1] 

The  bounds  calculated  in  the  previous  section  is  for  a  fixed  r  in  (0, 1],  i.e. ,  k(t) 
and  k(t),  where  the  argument  indicates  that  the  bounds  are  a  function  of  r.  To 
find  the  global  bounds,  the  algorithms  in  the  previous  must  be  performed  in 
r-space. 

Algorithm  2.2.9  Global  bounds 

1.  Set  Nr,  the  number  of  samples  in  (0, 1). 

2.  Generate  Nr  uniformly  distributed  random  numbers  in  (0, 1). 

3.  Initialise  the  set,  R  =  {ro,  ri,  r2, . . . ,  rjvr,  rjvr+i},  where  ro  =  0,  r*  for 
i  =  1,2, ,  Nr,  are  the  random  numbers  generated  in  #2,  and  r^r+ i  =  1. 

4-  For  all  rs  £  R  perform  the  following: 

•  Run  Algorithms  2.2.5,  and  2.2.6  for  <f>  or  2.2. 7  for  (f. 

•  Using  n(rs)  and  R(rs),  calculate  the  bounds  for  cf>k ,  i.e.,  ^max,  4>m 
or  f)min, 

5.  For  Vi  in  R,  compare  the  bounds  difference  with  either  r,_i  or  rj+i  in  R, 
generate  new  ri  towards  the  direction  where  the  bounds  improve. 

6.  If  |rjvr+i  —  ro|  is  less  than  a  tolerance,  then  stop.  Otherwise,  go  to  # 4 
with  the  updated  R. 

It  is  worth  to  note  that  the  suggested  algorithm  is  different  from  the  blind 
Monte-Carlo  random  search.  The  search  performs  through  random  samples  but 
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each  iteration  it  converges  towards  a  solution.  As  a  result,  the  search  space 
reduces  every  iteration. 

Remark  2.2.10  The  total  computational  cost  depends  on  how  many  times  the 
algorithms  are  executed  for  different  r  in  (0, 1] .  If  the  upper  bounds  could  be 
provided  with  less  computation,  the  proposed  algorithm  could  be  used  only  for 
finding  lower  bounds. 

The  proposed  four  algorithms  in  the  above  are  embarrassingly  parallel ,  i.e. 
all  sampling  evaluations  are  independent  from  each  other  and  it  can  be  easily 
implemented  on  parallel  computer  architectures.  Therefore,  it  is  a  perfect  prob¬ 
lem  to  be  solved  on  multi-core,  distributed  parallel  computer  nodes  and  GPU. 
For  some  examples  presented  in  the  next  chapter,  the  sampling  evaluation  part 
of  the  algorithms  is  running  on  NVIDIA  Tesla  C2050  GPU,  which  has  449  cores 
and  the  maximum  number  of  threads  per  block  is  1024.  The  algorithms  on  the 
GPU  is  implemented  using  CUDA-GPU  [27]. 

2.3  Probabilistic  Properties  of  the  Algorithms 

A  probabilistic  guarantee  of  the  performance  for  the  algorithm  finding  tj>  is  to 
be  derived  using  the  Chernoff  bound  [28] .  Similar  derivation  can  be  obtained 
for  the  other  parts  of  the  algorithms.  If  the  number  of  samples,  N,  satisfies  the 
following  inequality: 

(2.26) 

then  the  following  is  true: 

Probability  { \p(r)  -  p(r)  \  <  ec}  >  (1  -  Sc)  (2.27) 

for  ec  €  (0,  1)  and  5C  €  (0,  1),  where 

P(r)  =  (2.28) 
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Vi  is  the  connected  volume,  of  which  volume  size  is  assumed  to  be  strictly  greater 
than  zero,  where  all  values  of  </>fc  in  the  volume  are  greater  than  the  maximum 
lower  bound  found  by  the  algorithm,  Vt  is  the  total  volume  of  the  search  space, 


p(r)  = 


rk 

N’ 


(2.29) 


and  ni  is  the  number  of  samples,  which  belongs  to  Vi  .  Based  on  these,  we  can 
prove  the  following  theorem: 

Theorem  2.3.1  If  the  number  of  samples  is  equal  to  the  Chernoff  bound,  (2.26) 
and 

pW  £  Tv  +  {2 ^losf  •  (2'30) 

then  the  probability  of  finding  the  larger  lower  bound  than  the  current  one  using 
the  algorithm  is  greater  than  (1  —  (5C). 

Proof.  Let  N  equal  to  the  Chernoff  bound  in  (2.26), 


ec  = 


(2.31) 


From  (2.27), 


Probability  {max  [0,  p(r)  —  ec]  <  p(r)  <  min  [1,  p(r)  +  ec}}  >  1  —  5C  (2.32) 

Substituting  (2.29)  and  (2.31)  into  (2.32),  then  the  probability  of  the  following 
inequality  hold  is  greater  than  1  —  6C: 


max 


".pM-i/Ttogl 


<  ^  <  min[l.p(r)  +  ec]  (2.33) 


By  substituting  the  lower  bound  for  p(r),  (2.30),  into  (2.33), 


Probability  | N  max  ^0,  —  )  <  rq  [>  =  Probability  { 1  <  n{\  >  1  —  5C  (2.34) 
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where  Probability!  1  <  n{\  is  the  probability  of  hireling  at  least  one  sample, 
whose  4>k  is  greater  than  the  current  lower  bound.  ^ 

Note  that  if  the  true  probability  is  in  the  range: 

0<pW<^  +  y^1og^,  (2.35) 

then  the  bound  in  (2.34)  becomes 

Probability  {0  <  n;}  >  1  —  Sc  (2.36) 

i.e.  the  lower  bound  is  changed  from  1  to  0  in  the  left  hand  side  of  bracket.  As  a 
result,  in  this  case  we  cannot  say  anything  about  the  probability  for  improving 
the  lower  bound. 
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Chapter  3 


Results  and  Discussion 


3.1  “ Hello  World ”  Example 


Consider  the  following  <pk,  which  is  called  Rosenbrock  function: 


4>k{p)  =  100  (p2  -pf)2  +  (1  -Pif,  (3.1) 

where  it  is  assumed  that  there  is  no  uncertainty  in  the  initial  state,  and 


—7 r  <  Pi  <  n\/2 


(3.2) 


for  i  =  1  and  2.  Normalise  pi  such  that 


Pi  =  Pa  +  WpiSpi, 


(3.3) 


where 


r-s/2  —  - 


Pa  =  + 


tV2- 


(3.4a) 

(3.4b) 
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Figure  3.1:  Rosenbrock  function 


and 


Wp  = 


wPl 

0 


0 

wP2 


(3.5) 


The  function,  (3.1),  has  one  global  minimum  at  pi  =  1  and  P2  =  1,  whose 
function  value  is  zero,  and  one  global  maximum  at  p±  =  7ry/2  and  P2  =  —  7r, 
whose  function  value  is  about  52,365. 

For  Algorithm  2.2.5,  Pre-K  estimation,  N  =  1,000,  e  =  1  x  10“ 12 ,  E  = 
1  x  1012,  e  =  0.001,  and  the  calculated  kp  is  equal  to  1/52364.  For  Algorithms 
2.2.6  and  2.2.7,  Sf  =  1  +  10-6,  s  =  1.2,  e  and  E  are  the  same  as  ones  in  Pre-K 
estimation.  Finally,  these  have  to  be  run  for  different  r  in  Algorithm  2.2.9, 
where  Ay  =  21  and  the  tolerance  for  r,v.r+i  —  ?~o|  is  1  x  10-6.  The  MATLAB 
programmes  can  be  found  in  Appendix  A. 

From  Algorithm  2.2.6,  </^ax  ~  52,364.9  and  ^  ~  52,365.0,  where  the 
true  maximum  is  about  52,365.  On  the  other  hand,  from  Algorithm  2.2.7, 
<bk  .  «  —0.00001  and  6^-  «  0.063,  where  the  true  minimum  is  equal  to  0. 
Because  of  the  probabilistic  nature  of  the  upper  bounds,  they  cannot  be  trusted 
with  100%  confidence.  However,  the  lower  bounds  are  extremely  close  to  the 
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Table  3.1:  Bounds  for  the  minimum  of  high-dimensional  Rosenbrock  function, 
where  the  true  minimum,  0m;n,  is  equal  to  0. 


n 

2 

3 

4 

5 

6 

0min 

0.063 

0.0906 

0.1877 

0.2513 

-0.0664 

0  • 
—min 

-0.00001 

-0.00046 

-0.00845 

-0.0056 

-0.2539 

true  values  in  this  case. 

High-dimensional  Rosenbrock  function  can  be  defined  as 

n—  1 

4>k{p)  =  100  (Pi+i  ~  PiY  +  (!  -Pi)2  (3-6) 

i= 1 

where  n  is  the  dimension  of  the  uncertain  space.  The  global  minimum  occurs  at 
all  Pi  =  1  for  i  =  1,2, ...  ,n  and  the  minimum  is  zero.  For  different  n,  the  bounds 
for  the  minimum  are  shown  in  Table  3.1.  The  gap  between  the  upper  and  the 
lower  bounds  becomes  larger  as  the  dimension  of  the  uncertain  parameter  space 
being  larger.  It  eventually  fails  at  n  =  6,  where  the  upper  bound  is  smaller  than 
0.  To  resolve  this,  the  number  of  samples,  N,  is  increased  from  1000  to  2000. 
And,  the  bounds  obtained  is  </>min  =  0.187  and  6  .  =  —0.0683. 

The  relation  between  the  bounds  gap  and  the  specified  N  is  not  known  in 
general.  That  is,  how  N  must  be  chosen  in  order  to  make  sure  that  the  bounds 
is  correct  is  unknown.  In  addition,  the  required  N  increases  exponentially  as 
the  search  space  dimension  increases.  For  example,  from  (2.27),  N  equal  to 
2000  gives  99.99%  confidence  in  finding  p(r)  greater  than  0.05%.  However,  p(r) 
around  the  global  minimum  reduces  quickly  as  the  search  dimension  increases. 
Practically,  up  to  a  certain  range  of  n,  it  can  be  effectively  overcome  with 
increased  samples  and  by  parallelising  the  algorithms  over  distributed/GPU 
processing  units. 
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Figure  3.2:  Bounds  for  x\^ 


3.2  Oscillatory  state 


The  following  example  is  a  discretised  rotating  system  [20]: 


1  p 

Xi,k 

%2,k+l 

r- - ^ 

V1+P 

-p  1 

%2,k 

(3.7) 


In  the  original  algorithm  in  [20],  the  uncertainty  must  be  in  a  polynomial  for¬ 
mat  and  a/1  +pl  was  used  instead  of  a/1  +  p2.  Here,  it  does  not  need  to  be 
polynomial  and  the  original  form,  a/1  +  p2 ,  is  used.  The  intervals  for  the  initial 
state  and  the  uncertain  parameters,  p,  are  given  as  follows:  0.9  <  £o,i  <  1.1, 
0.9  <  £0,2  <  1-1,  and  0.45  <  p  <  0.55. 

Algorithm  2.2.5  firstly  calculates  the  pre  upper  bound  for  ]</>],  1/k,  and  set 
s  =  2.  Secondly,  the  bounds  for  max(i))  and  min(<)>)  are  obtained  by  Algorithms 
2.2.7  and  2.2.6.  For  many  practical  systems,  extreme  (both  maximum  and  the 
minimum)  likely  occurs  at  the  boundary  of  at  least  one  uncertain  parameters. 
Hence,  frequently,  sweeping  algorithm  over  r,  Algorithm  2.2.9,  is  not  required 
but  the  running  algorithms  for  r  =  1  is  sufficient. 

Figure  3.2  shows  the  bounds  for  x±}k,  where  k  =  1,  2,  . . . ,  39,  40.  The  upper 
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and  lower  bounds  of  the  maximum  and  the  minimum  for  r  =  1  shown  in  Figure 
3.2  are  very  close  to  each  other.  All  trajectories  from  random  simulations,  where 
the  samples  are  from  all  r  in  (0, 1],  are  well  bounded  by  the  estimated  bounds. 

The  calculation  time  becomes  longer  as  k  increases.  A  simple  solution  is 
resetting  the  initial  starting  step  once  in  a  while  as  presented  in  [20].  This, 
however,  will  reduce  the  computational  time  with  the  cost  of  less  tight  bounds. 
On  the  other  hand,  the  state  bounds  prediction  in  practice  may  not  require  very 
large  k. 


3.3  ErbB  Signalling  Pathways 

Epidermal  growth  factor  receptor  (ErbB)  related  pathways  are  among  the  most 
extensively  studied  biological  signalling  networks  [14].  Abnormality  of  ErbB 
signalling  pathways  cause  various  human  cancers  [29,  30,  31].  In  [14],  an  ErbB 
mathematical  model  including  13  known  ErbB  ligands,  Epidermal  Growth  Fac¬ 
tor  (EGF)  and  Heregulin  (HGF)  and  Erk  and  Akt  pathways  are  presented.  It 
has  504  states,  828  reactions  and  226  kinetic  parameters.  The  set  of  504  differ¬ 
ential  equations  is  extracted  from  the  simbiology  model  [14].  It  is  known  that 
this  model  is  only  valid  up  to  a  few  hours  and  it  is  not  necessary  for  this  system 
to  be  stable  for  infinite  time  interval.  As  long  as  the  states  remain  in  a  certain 
bound,  the  network  works  perfectly  as  it  should  be.  Hence,  the  required  robust¬ 
ness  analysis  is  obtaining  the  information  about  how  the  future  state  bounds 
propagate  with  respect  to  the  uncertain  parameters. 

One  of  the  interesting  biological  features  found  in  [14]  using  the  model  is 
that  parametric  sensitivities  of  the  dynamical  system  strongly  depend  on  in¬ 
put  condition.  This  could  be  the  reason  that  it  provides  so  diverse  responses. 
Parametric  uncertainties  are  introduced  for  those  226  kinetic  parameters.  The 
uncertainty  ranges  are  set  to  ±10%  from  the  nominal  values.  Among  the  several 
input  conditions,  the  robustness  is  tested  for  the  case  of  EGF  equal  to  5nM.  The 
bounds  for  phospholilated  ErbBl,  i.e.  p-ErbBl,  is  shown  in  Figure  3.3.  Again, 
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Figure  3.3:  Bounds  for  p-ErbBl  with  respect  to  uncertainties  in  226  kinetic 
parameters 

the  upper  and  lower  bounds  for  the  maximum  and  the  minimum  are  very  close 
to  each  other  and  only  the  upper  bounds  for  both  are  indicated.  All  trajectories 
from  random  simulations  are  well  bounded  by  the  estimated  bounds.  For  the 
bounds  calculation,  r  is  fixed  to  1  and  for  the  random  simulations,  r  is  randomly 
selected  between  0  and  1. 

As  the  size  of  the  system  is  relatively  large,  the  integration  part  is  imple¬ 
mented  over  CUDA-GPU  [27]  and  the  programmes  are  shown  in  Appendix  B. 

3.4  Inverted  Pendulum:  Hybrid  System 

A  switching  controller  for  inverted  pendulum  stabilisation  is  shown  in  [32].  A 
simplified  version  of  the  system  is  given  by 

9  =  p  sin  9  —  u  cos  $,  (3.8) 

where  8  is  the  angle  of  pendulum  measured  from  the  upright  position,  p  is  the 
uncertainty  caused  by  some  physical  parameters,  and  u  is  the  control  input. 
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Figure  3.4:  Controller  switching  phase  examples:  0  for  waiting;  1  for  the  energy 
dissipation  control;  and  2  for  feedback  linearisation  control. 


The  total  energy,  E,  including  kinetic  and  potential  energy  is  given  by 

E  =  +  (cos  9  -  1) 


The  controller  proposed  in  [32]  has  the  following  switching  behaviours: 

•  Energy  dissipation  phase  (phase  1):  if  \E\  >  e,  where  e  is  a  positive  real 
number,  then 


sign(i£)0 

l  +  l^l 


(3.9) 


•  Waiting  phase  (phase  0):  if  \E\  <  e  and  \6\  +  \6\  >  6, 


u  =  0 


(3.10) 


Feedback  linearisation  control  phase  (phase  2):  if  \E\  <  e  and  \0\  +  \6\  <  S , 


29  +  9  +  sin  6 
cos  9 


(3.11) 


The  ranges  for  the  initial  values  are  set  to:  |#(0)|  <  19°,  |0(O)|  <  20°/s, 
the  uncertainty  range  is  given  by  0.1  <  p  <  1.9,  i.e.  ±90%  uncertainty  from 
the  nominal  value,  1,  and  e  and  5  are  set  to  0.1  and  0.8,  respectively.  An 
example  of  the  controller  phase  switching  history  is  shown  in  Figure  3.4.  The 
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Figure  3.5:  Bounds  for  6{tk) 


system  dynamics  is  highly  nonlinear  because  of  its  inherent  nonlinearities  and 
the  switching  control. 

The  estimated  bounds  obtained  with  r  =  1  are  shown  in  Figure  3.5.  Al¬ 
though  it  only  requires  to  calculate  the  bounds  for  one  of  either  sign  as  the 
system  is  symmetric,  for  the  demonstration  purpose,  both  bounds  are  calcu¬ 
lated.  At  t  =  10s,  the  bounds  are  between  ±140°  but  most  of  the  trajectories 
found  by  Monte-Carlo  simulations  converge  to  zero.  In  fact,  the  Monte-Carlo 
simulation  method  finds  only  one  trajectory  but  still  far  from  the  bound  at 
t  =  10s. 

This  clearly  demonstrates  the  advantage  of  the  proposed  bound  algorithm 
over  the  blind  Monte-Carlo  simulations.  The  number  of  samples,  N,  for  this 
example  is  set  to  1024x10  and  the  one  for  Monte-Carlo  simulations  is  twice  more 
than  N  used.  The  calculation  time  of  the  presented  algorithm  for  each  instance 
is  less  than  0.5s.  Monte-Carlo  simulations  takes  significantly  longer  time,  about 
3  minutes,  which  would  vary  depending  on  the  number  of  samples.  The  Monte- 
Carlo  random  simulation  cannot  find  any  solution  closer  to  the  lower  bounds  at 
t  =  10s. 
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Figure  3.6:  The  experiment  was  performed  in  the  rectangular  box.  The  arrows 
are  the  average  current  velocity  at  each  location  for  seven  different  dates  from 
11th  to  17th  March  of  2011  [33].  The  red  circles  in  the  inset  box  are  the 
position  information  from  GPS  and  the  black  thick  line  is  the  path  obtained  by 
the  navigation  system. 

3.5  Underwater  Glider  Navigation  System 

In  this  example,  Kinematics  for  underwater  glider  is  derived,  and  the  bounds 
for  the  future  position  of  underwater  glider  are  to  be  calculated  using  the  er¬ 
ror  model  including  sea  current  and  uncertainties  in  physical  parameters.  All 
experimental  data  is  obtained  by  Korea  Institute  of  Ocean  Science  &  Technol¬ 
ogy  (KIOST)  in  the  area  indicated  in  Figure  3.6. 

3.5.1  Kinematics 

The  relation  among  the  position  vector  of  the  underwater  glider,  r^,  the  sensor 
position  in  the  body-coordinates  (£>),  rs,  and  the  sensor  position  in  the  reference- 
coordinates  (7 Z),  r,  is  shown  in  Figure  3.7  and 


(3.12) 


where  is  the  position  vector  expressed  by  two  vectors  in  the  right  hand  side, 
and  B  in  the  superscripts  is  indicated  that  both  vectors  are  expressed  in  the 
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Figure  3.7:  The  origins  of  the  reference-coordinates  and  the  body-coordinates 
are  indicated  by  O  and  B ,  respectively.  The  body  coordinates  is  attached  to 
the  glider  and  rs  is  the  sensor  position  vector  relative  to  the  origin  of  the  body- 
coordinates.  g  is  the  gravitational  acceleration  vector. 


body-coordinates.  Take  the  time-derivative  of  (3.12)  and  obtain  the  velocity, 


vb  = 


—  us  x  r 


B 

S  > 


(3.13) 


where  rf  is  zero  because  the  glider  is  assumed  to  be  a  rigid-body,  and  the 
following  transport  theorem  is  used  in  the  derivative  [34].  Take  one  more  time- 
derivative  of  (3.13) 


ab=dtVb=dfir~UJXT‘ 


—  uj  x  (u>  x  ra) , 


(3.14) 


where  the  superscript,  B ,  in  the  terms  in  the  right  hand  side  is  dropped  for 
brevity,  and  they  are  all  expressed  in  the  body-coordinates. 


The  following  quantity,  a”cc,  is  the  measurement  of  the  3-axis  accelerometer 
in  the  glider: 

aL  =  +  cn  (3-15) 

where  9 ,  <j>)  is  the  direction  cosine  matrix,  which  changes  the  coordinates 

from  the  reference  to  the  body  with  the  rotating  sequence  given  by  yaw(^), 
pitch(0),  and  roll(^),  g7^  =  [0,  0,  —  g]T ,  the  gravitational  acceleration  expressed 


31 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


in  the  reference  coordinates,  and  g  =  9.821  m/s2.  Substitute  d2r/dt2  in  (3.14) 
by  the  expression  in  (3.15) 

afc  =  =  afcc  -  C%gK  -  u>  x  rs  -  u>  x  (u>  x  rs) .  (3.16) 

Rotational  angular  velocity  of  underwater  vehicle  is,  in  general,  very  slow,  i.e. 

1 1 n;||  and  ||d>||  are  much  smaller  than  ||g||.  In  the  experiment  the  magnitudes  of 
two  terms  are  at  least  100-times  smaller  than  the  magnitude  of  the  gravitational 
acceleration.  Hence, 


ab  «  afcc  -  C%gK.  (3.17) 

After  transforming  the  coordinates  of  (3.17)  into  the  reference-coordinates,  inte¬ 
grating  it  with  respect  to  time  provides  the  velocity  and  the  position  expressed 
in  the  reference  frame. 

3.5.2  Error  Model 

Exocetus  Coastal  glider  has  a  built-in  navigation  algorithm  and  it  records  vari¬ 
ous  navigation  data.  The  following  derivations  are  based  on  the  log-file  obtained 
from  the  glider  navigation  system.  The  3- Axis  Acceleration  measurements  ex¬ 
pressed  in  the  body-coordinates,  afcc,  are  directly  available  from  the  log-file  and 
it  can  be  expressed  as  follows: 

aacc  =  r  +  C tj(^>,  9 ,  <f>) gK  +  bacc  +  nacc,  (3.18) 

where  r  is  the  measured  quantity  by  the  calibrated  accelerometer  if  there  is  no 
stochastic  noise  and  bias  error,  bacc  is  the  accelerometer  sensor  bias  error,  nacc 
is  the  white  noise,  and  the  accelerations  caused  by  the  control  forces  and  the 
external  disturbances  are  all  included  in  r. 

Usually,  in  the  estimation  problem,  r  is  assumed  to  equal  to  r,  which  pre¬ 
sumes  infinite  resolution  and  infinitesimal  time  constant,  i.e.  a  perfect  sensitiv- 
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ity.  In  reality,  including  the  limitations  in  the  resolution  and  the  finite  response 
time,  and  some  unknown  error  sources,  r  will  be  different  from  r  [35]. 

Attitude  Error:  The  navigation  algorithm  provides  the  attitude  information, 
which  includes  some  errors,  Sip,  69,  and  Sep,  then 

C\  =  C%, (Sip,  80,54 )  C%(4,0A)i  (3.19) 

where  B'  is  the  body-coordinates  that  the  navigation  system  calculated,  ip,  9,  <p 
are  roll,  pitch  and  yaw  angles  returned  from  the  navigation  algorithm,  and 
Sip,  69,  6cp  are  the  attitude  angles  of  the  true  body-coordinates,  B,  with  respect 
to  the  calculated  body-coordinates,  B' . 

Accelerometer  Resolution  Limit:  The  acceleration  of  the  glider,  (3.17),  is 
calculated  as  follows: 


ab  =  afcc  -  (3.20) 

Once  the  bias  correction  is  completed,  the  position  information  can  be  obtained 
by  integrating  the  corrected  accelerometer  output  twice.  The  calculated  path  is 
shown  in  the  inset  of  Figure  3.6  indicated  in  the  black  solid  line.  The  starting 
position  was  initialised  by  the  Grobal  Positioning  System  (GPS)  and  it  can  be 
considered  as  an  accurate  starting  position  initialisation.  However,  because  of 
acceleration  below  the  sensor  resolution  is  accumulated  during  the  diving,  the 
position  at  the  end,  just  before  updated  by  GPS  again,  has  several  km  error 
after  about  3  hours  operation  under  the  water.  The  main  cause  of  this  error 
might  be  the  current  as  the  final  position  is  almost  exactly  drifted  towards  the 
directions  where  the  average  current  flows.  The  similar  drifts  are  observed  in 
another  three  divings  [36]. 

Sea  Current  Uncertainty:  Av  is  the  relative  velocity  of  the  glider  in  the 
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horizontal  plane,  i.e. 


T;crt 

vx 

Av  = 

C,crt 

% 

- 

Vy 

0 

0 

(3.21) 


v and  u£jrt  are  the  estimated  sea  current  velocity  in  the  east  or  the  north  direc¬ 
tion,  respectively,  whose  magnitude  is  around  lkt  («  0.5  m/s)  in  the  maximum, 


<=«Ert(l  +  ^E),  (3.22a) 

<=^‘(1  +  ^),  (3.22b) 

Ujjf*  and  Uj/rt  are  the  true  sea  current  velocity  in  the  east  or  the  north  direction, 
respectively,  and  Sve  and  5i>n  are  the  uncertainties. 

Parameters  Uncertainty  in  Drag:  To  quantify  how  much  acceleration  would 
be  generated  from  the  current,  the  following  equation  is  used  [37]: 

D=  X-p  ||Av||2  AdCd(1  +  6D),  (3.23) 

where  p  is  the  sea  water  density,  which  is  around  1020  kg/m3,  Ad  is  the  cross 
section  area  of  the  glider  hull,  which  is  about  824  cm2,  and  Cd  is  the  drag 
coefficient  set  to  0.4,  which  is  adopted  from  [37].  Unlike  the  example  shown 
in  [37],  the  glider  has  a  small  wing  and  the  induced  drag  is  ignored.  As  these 
physical  parameters  are  estimated  using  a  similar  size  glider,  the  calculated 
drag  has  some  error  and  it  is  represented  by  SD.  With  the  nominal  values,  i.e. 
SD  =  0,  the  drag  is  about  4N  when  vx  =  0,  and  vy  =  0.  It  corresponds  to  the 
acceleration  of  0.04  m/s2  («  4N/109kg).  The  acceleration,  0.04  m/s2,  is  the 
possible  maximum  value  when  the  glider  is  stationary.  The  acceleration  caused 
by  the  current  would  be  much  smaller  than  the  maximum  as  the  glider  usually 


34 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


does  not  fly  in  the  exact  opposite  to  the  current.  Then,  these  tiny  current  effect 
would  be  completely  blocked  by  the  sensor  noise  and  the  external  fluctuations. 


3.5.3  Longitude/Latitude  Bound  Estimation 

The  acceleration  caused  by  the  drag  with  considering  the  sea  current  is  given 

by 

w  «  —Cjg,  Ar  =  — eD,  (3.24) 


where  m  is  the  mass  of  the  glider  and  &d  is  the  unit  vector  towards  the  cur¬ 
rent  direction.  Sea  current  near  the  surface  is  mainly  caused  by  wind  and  the 
temperature  differences.  The  vertical  direction  is  negligible  in  the  near  surface 
depth  compared  to  the  magnitudes  of  the  horizontal  direction  components  [38] . 
Hence, 


Av 


= 


||Av| 

0.  otherwise 


,  for  ||  Av ||  >  0 


(3.25) 


Scenario  # 1 :  Set  the  uncertain  parameters  as  follows: 

-5°  < 6 if;  <  5°, 

(3.26a) 

-5°  <66  <  5°, 

(3.26b) 

-5°  <6(j><  5°, 

(3.26c) 

—0.2  [m/s]  <5v e  <  0.2  [m/s], 

(3.26d) 

—0.2  [m/s]  <5i>n  <  0.2  [m/s], 

(3.26e) 

-0.1  <5D  <  0.1, 

(3.26f) 

and  the  nominal  values  of  the  sea  current  velocity  are:  Dgrt  = 

—0.1165  m/s  and 
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Figure  3.8:  Paths  after  1  minutes  by  Monte  Carlo  simulations  are  indicated  by 
green  lines  and  the  lower  bounds  are  indicated  by  dashed  lines.  The  dot  at  the 
end  of  each  path  is  the  final  location  of  the  glider.  All  dots  are  inside  the  box 
defined  by  the  dashed  line. 

rV  =  0.3948  m/s.  The  bounds  for  the  possible  longitude  and  latitude  coordi¬ 
nates  after  60s  are  calculated  with  N  =  100  and  r  is  fixed  to  1.  The  four  lower 
bounds  are  indicated  in  Figure  3.8.  All  paths  start  at  the  same  coordinates 
and  the  final  position  for  each  random  simulation  is  indicated  by  blue  dot.  All 
final  positions  are  well  bounded  by  the  lower  bound.  The  upper  bounds  are  not 
indicated. 

Scenario  As  the  uncertain  ranges  in  Scenario  #1  are  too  big  for  longer 
time  horizon,  i.e.,  the  paths  will  spread  broad  region  with  longer  final  time, 
reduce  the  uncertain  parameter  ranges  by  10  as  follows: 


-0.5°  <5i)  <  0.5°, 

(3.27a) 

-0.5°  <59  <  0.5°, 

(3.27b) 

-0.5°  <5(p  <  0.5°, 

(3.27c) 

0.02  [m/s]  <<5ue  <  0.02  [m/s], 

(3.27d) 

0.02  [m/s]  <5v n  <  0.02  [m/s], 

(3.27e) 

-0.01  <5D  <  0.01. 

(3.27f) 
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Figure  3.9:  Paths  after  1  hour  by  Monte  Carlo  simulations  are  indicated  by 
green  lines  and  the  lower  bounds  are  indicated  by  dashed  lines.  The  dot  at  the 
end  of  each  path  is  the  final  location  of  the  glider.  All  dots  are  inside  the  box 
defined  by  the  dashed  line. 


The  final  time  is  set  to  1  hour.  The  lower  bounds  shown  in  Figure  3.9  confine 
tightly  all  the  possible  positions  after  1  hour.  The  number  of  Monte  Carlo 
simulation  is  500. 
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Chapter  4 


Conclusions 


Algorithms  for  calculating  state  bounds  for  general  nonlinear  systems  with  un¬ 
certain  initial  states  and  parameters  are  developed.  The  algorithms  combine 
/x-formulation  for  optimisation  problem  and  pseudo-LFT  for  LFT-free  formula¬ 
tion.  The  algorithms  have  several  advantages  including:  1)  no  effort  to  obtain 
an  LFT  form;  2)  easy  to  parallelise  on  distributed  computers;  and  3)  applicable 
to  many  types  of  functions  including  the  one  having  finite  number  of  disconti¬ 
nuity.  The  algorithms  are  applied  to  multi-dimensional  Rosenbrock  function;  a 
simple  oscillatory  nonlinear  discrete  system;  a  high- dimensional  biological  model 
for  ErbB  signalling  pathways;  a  hybrid  system  with  discontinuity;  and  finally, 
navigation  error  propagation  of  underwater  glider. 

It  is  highly  desirable  to  have  numerically  efficient  algorithms  to  estimate  state 
bounds  for  general  nonlinear  systems.  Especially,  the  input-output  robustness 
analysis  with  respect  to  various  parametric  perturbations  are  one  of  the  main 
interest  in  the  robustness  of  biological  networks;  autonomous  mobile  robots 
operating  in  uncertain  environment  requires  to  predict  the  future  state  bounds 
in  order  to  plan  or  re-plan  its  behaviour;  predicting  a  group  of  space  debris  is 
very  important  for  the  safety  of  space  mission.  The  suggested  algorithms  could 
be  the  main  tool  to  analyse  such  complex  nonlinear  systems  with  uncertainties. 

As  demonstrated  by  the  multi-dimensional  Rosenbrock  functions,  however, 
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the  gap  between  the  lower  bounds  and  the  upper  bounds  increases,  in  general, 
as  the  search  dimension  increases.  This  is  the  inherent  complexity  of  NP-hard 
problem,  which  cannot  be  overcome  unless  NP=P.  On  the  other  hand,  up  to 
certain  dimensions,  tighter  gap  can  be  obtained  with  much  shorter  computation 
time  by  parallelising  the  algorithms.  We  demonstrated  the  possibility  of  parallel 
processing  in  real-time  using  GPU.  In  future,  it  will  be  possible  to  perform  real¬ 
time  optimisation  for  predicting  uncertain  bounds  and  adjusting  control  gains 
to  reduce  the  bounds  with  fast  massive  parallel  processing  computation. 

Finally,  it  is  worth  to  note  that  a  similar  algorithm  to  the  proposed  ones  was 
presented  in  [39],  which  is  called  Luus-Jaakola  algorithm,  and  its  convergence 
property  was  shown  in  [40] .  It  uses  random  sampling  and  re-sanrpling  based  on 
the  values  from  the  random  samples.  It  was  presented  in  1973  and  its  usages 
were  limited  to  small  size  problems  because  of  a  lot  restricted  computational 
power  in  the  past.  The  re-sampling  method  could  be  adopted  in  the  proposed 
algorithm  to  reduce  the  computational  cost. 
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Appendix  A 

Hello  World  example 

A.l  main_get_max_all.m 


2  %  Main:  States  bound  calculation  for  many  k  steps 

3  % - 

4  % 

5  %  See  #1,  #2  and  #3  below 

6 

7  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

8  %  University  of  Glasgow 

9  %  Biomedical  Engineering 

10  %  RobustLab:  http://www.robustlab.org 

n 

12  %  Change  Log: 

13  %  1.  15th  January  2014 

14  %  To  solve  optimisation  problem 

15  %  2.  15th  Feburuary  2014 

16  %  To  correct  algorithm  for  finding  solutions  inside  boundary 

17  %  3.  1st  April  2014 

is  %  To  find  the  maximum  or  the  minimum;  remove  the  positiveness  condition 

19 

20  clear  all; 

21 

22  %% 

23  % - 

24  %  1 .  Algorithm  Parameters  that  you  may  want  to  adjust 

25  % - 

26 

27  %  number  of  samples  in  one  Period 

28  num.steps  =  1; 

29  num_delta_sample  =  1000; 

30  %  lower  with  less  computation  but  higher  risk  to  be  failed 

31  %  higher  with  less  risk  but  more  computational  burden 

32 

33  saf  ety_f  actor_max  =  1  +  le— 6; 

34  saf  ety_f  actor_min  =  1  +  le— 6; 

35  %  always  greater  than  or  equal  to  1 

36  %  smaller  would  give  tighter  upper  bounds  with  higher  risk  to  be  failed, 

37  %  larger  would  have  less  possbility  to  be  failed  but  less  tight  upper  bounds 

38 

39  %  find  minimum  or  maximum 
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40  fincLmin  =  true;  %  true  for  minimum;  false  for  maximum 

41 

42  %  Range  of  k,  where  k  for  k*delta_c 

43  box.init  =  [le— 12  lel2];%[— 0.5  l]*le3; 

44  %  For  box_init  =  [le— 12  lel2] ; 

45  %  you  should  know  roughy  where  the  k.optimal  would  be. 

46  %  i.e,  here  I  am  quite  sure  that  le— 2  <  f  <  lel2 

47  %  if  your  calculation  failed,  then  you  may  adjust  this 

48 

49  %  Use  adaptive  box  &  samplings 

50  use.adapt  ive.sampling  =  false; 

51  num_replace  =  f loor (num_delta.sample*0 . 8+0 . 5) ; 

52 

53  %  - 

54  %  2.  the  function  f(x,p)  related  setting 

55  %  - 

56  %  parameter  bounds 

57  x_dim  =  2; 

58  xL  =  —  pi*ones  ( x_dim,  1 )  ; 

59  xU  =  sqrt (2) *pi*ones (x.dim, 1) ; 

60 

61  s.scale  =  1.2; 

62 

63  %  uncertain  parameters  dimension 

64  state.dim  =  length  (xL); 

65 

66  w_x  =  (xU— xL)  /2; 

67  x_centre  =  xL+w_x; 

68 

69  w.x.prev  =  w_x; 

70  x.centre.prev  =  x.centre; 

71 

72  w_x.org  =  w_x; 

73  x_centre.org  =  x.centre; 

74 

75  %% 

76  %================================================== 

77  %  No  need  to  change  from  here 

78  %================================================== 

79  %  find  maximum  of  |f| 
so  tic 

si  c_max  =  0; 

82  [x.upper.max ,  x.lower.max,  sol.opt,  weight.opt]  =  ... 

83  mult  i_cal_state_bounds  (  [  ]  ,  ... 

84  box.init ,  num.steps  ,  num.delt a.sample ,[],  state.dim,  ... 

85  [],[],[],  x.centre ,  w.x ,  [  ]  ,  ... 

86  [],[],  num.replace ,  [  ]  ,  c.max ,  find_min)  ; 

87  opt.x  =  x.centre  +  weight.opt  *sol_opt  (:).  *w_x  (:)  ; 

88  fprintf (' Optimal  Sol:  x  =  ('); 

89  for  idx=l :  state.dim 

90  fprintf  (  1  %4 . 2f  '  ,  opt.x  (idx)  )  ; 

91  end 

92  fprintf  (  '  )  \n  '  )  ; 

93 

94  x.opt  =  x.centre  (:)'  +  ( sol.opt  (:).  *w_x  (:))'  ^weight.opt ; 

95  opt.cost  =  get.f  .delta  ( sol.opt ,  x.centre ,  w_x*weight_opt ,  0 )  ; 

96  fprintf (' Optimal  Cost  =  %10 . 5f\n ', opt.cost— c.max) ; 

97  fprintf (' Cost  bound  =  [%10.5f,  %10.5f]\n',  ... 

98  x_lower_max- c_max ,  x.uppe r .max— c.max )  ; 

99 

100  cal.time  =  toe; 

101  fprintf  (  ' - \n  '  )  ; 
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102  fprint f (' Calculation  Time  =  %6.3f  [s  ]  \n  1 ,  cal.time)  ; 

103  fprintf  (  ' - \n '  )  ; 

104 

105  f_abs_max  =  opt.cost; 

106 

107  %  find  maximum  of  |  f +f  _abs_max  *  2  | 

108  tic 

109  c_max  =  f.abs_max*s_scale; 

no  [x_upper_max,  x.lower.max,  sol.opt,  weight.opt]  =  ... 

m  mult  i_cal_state_bounds  (  [  ]  ,  ... 

112  box.init ,  num.steps  ,  num.delt a.sample ,  [  ]  ,  state.dim,  ... 

113  []/^  []/■  []/■  x_centre ,  w_x ,  [  ]  ,  ... 

114  [],  [],  num.replace ,  [  ]  ,  c_max ,  find_min)  ; 

115  opt_x  =  x.centre  +  weight_opt *sol_opt  ( : )  .  *w_x  ( : )  ; 

116  fprintf  (' Optimal  Sol:  x  =  ('); 

117  for  idx=l :  state_dim 

ns  fprintf (' %4 . 2f  '  ,  opt _x  ( idx)  )  ; 

119  end 

120  fprintf  (  '  )  \n  '  )  ; 

121 

122  opt_cost  =  get_f -delta  ( sol_opt ,  x_centre ,  w_x*weight_opt ,  0 )  ; 

123  fprintf  (' Optimal  Cost  =  %10 . 5f\n opt.cost ) ; 

124  fprint  f  (  '  Cost  bound  =  [%10.5f,  %10 . 5f  ]  \n  ' ,  x_lower_max— c_max,  ... 

125  saf ety_f act or_max*x_upper_max- c_max )  ; 

126 

127  cal.time  =  toe; 

128  fprintf  (  ' - \n  '  )  ; 

129  fprint f (' Calculation  Time  =  %6.3f  [s]  \n' ,  cal.time)  ; 

130  fprintf  (  1 - \n  '  )  ; 

131 

132  %  find  minimum  of  |f— c_max*2| 

133  tic 

134  c_max  =  -f_abs_max*s.scale; 

135  [  x_upper_max ,  x_lower_max,  sol.opt,  weight.opt]  =  multi_cal_state_bounds  (  [  ]  , 

136  box_init ,  num_steps  ,  num_delt  a.sample ,  [  ]  ,  state_dim,  ... 

137  [  ]  ,  [  ]  ,  [  ]  ,  x.centre ,  w_x ,  [  ]  ,  ... 

138  [],  [],  num_replace ,  [  ]  ,  c_max ,  find_min)  ; 

139  opt_x  =  x.centre  +  weight.opt  *sol_opt  (:).  *w_x  (:)  ; 

140  fprint  f  (' Optimal  Sol:  x  =  ('); 

141  for  idx=l :  state.dim 

142  fprintf  (  '  %4 . 2f  '  ,  opt.x  (idx)  )  ; 

143  end 

144  fprintf  (  '  )  \n  '  )  ; 

145 

146  opt.cost  =  get.f  .delta  ( sol.opt ,  x.centre ,  w_x*weight_opt ,  0 )  ; 

147  fprintf  (' Optimal  Cost  =  %10 . 5f\n  opt.cost )  ; 

148  fprint  f  (' Cost  bound  =  [%10.5f,  %10 . 5f  ]  \n  ' ,  x.lower.max+c.max,  ... 

149  saf ety.f act or_min*x_upper_max+c_max )  ; 

150 

151  cal.time  =  toe; 

152  fprintf  (  ' - \n  '  )  ; 

153  fprint f (' Calculation  Time  =  %6.3f  [s] \n cal.time) ; 

154  fprintf  (  ' - \n  '  )  ; 
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A. 2  multi_cal_state_bounds.m 
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function  [x_upper_max_global,  x.lower _max_global ,  extreme.delta,  extreme.w] 
mult  i_cal_state_bounds  (k.step,  box.init,  ... 

num.steps,  num.delta.sample ,  delta.dim,  state_dim,  which.x,  ... 
p.centre ,  w_p ,  x_cent re ,  w_x ,  saf ety.f actor ,  .  .  . 

f.real.prev,  delta.px.prev,  num.replace,  ... 
reset.steps,  c_max,  find_min) 

%  multiple  box  size 
Nr  =  21; 

w_box_size  =  linspace  ( 0 ,  1 ,  Nr)  ; 
max_iter  =  10000; 

min.iter  =  10; 

w.box.min.size  —  le— 6; 


w_x  =  w_x  (  : )  '  ; 

[zero.cost,  ]  =  get_f_delta  (zeros  (size  (w_x)  )  ,x_centre,  w_x,  c_max)  ; 


%%  calcualte  upper 

not  _c  on  verge  _sw 

x_upper_max_global 

x_lower_max_global 

extreme-delta 

extreme.w 


and  lower  bounds  for  the  maximum 
=  true; 

=  -inf; 

=  zero.cost; 

=  zeros (size  (w.x  ( : )  '  ) ) ; 

=  0; 


current.iter  =  1; 
prev.opt.x  =  zeros (size (w.x )) ; 
count.converge  =  0; 
w.prev.box  =  []; 


while  not.converge.sw 


sol_opt_all 
opt.cost.all 
opt_cost_all  (1) 


=  zeros  (length  (w.box.size )  ,  state.dim)  ; 
=  zeros ( length  (w.box.s i ze ), 1 ) ; 

=  zero.cost; 


current.w.range  =  w.box.size  (end)  —  w.box.size  (1 )  ; 


for  wdx  =  2 : length (w.box.si ze ) 

w.x.current  =  w_x*w_box_size  (wdx)  ; 

[x_upper_max,  x_lower_max,  sol.opt]  =  cal.state.bounds  (  [  ]  ,  ... 

box.init ,  num.steps ,  num_delta_sample ,  [],  state.dim,  ... 

[  ]  ,  [  ]  ,  [  ]  ,  x.centre ,  w.x.current ,  [  ]  ,  ... 

[]f  t]^  num_replace,  [],c.max,  find_min)  ; 

opt.cost  =  0; 
if  "isempty ( sol.opt ) 

[opt.cost,  ]  =  get.f.delta  (sol.opt ,  x.centre,  w.x.current ,  c.max)  ; 

end 


if  x.upper.max  >  x_upper_max_global 
x_upper_max_global  =  x_upper_max; 

end 

if  x.lower.max  >  x_lower_max_global 

x_lower_max_global  =  x_lower_max; 
extreme.delta  =  sol.opt  ( :  )  ' ; 
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59 


ext  reme.w 


=  w_box_size  (wdx)  ; 


60  end 

61 

62  sol_opt_all  (wdx,  : )  =  sol_opt(:)'; 

63  opt_cost_all  (wdx)  =  opt_cost; 

64 

65  end 

66  cost_dif  f  _sgn  =  sign  (dif  f  (opt_cost_all )  )  ; 

67  nbhd.idx  =  [2  :  length  (w_box_size )  ]  +  cost.dif  f  _sgn  ( : )  ' ; 

68  if  nbhd_idx  (end)>length  (w_box_size) 

69  nbhd.idx  (end)  =  length  (w.box.si ze )  ; 

70  end 

71 

72  figure  (1);  elf; 

73  [aa,bb]  =  sort  (w.box.si ze )  ; 

74  plot  (aa,  opt.cost.all  (bb)  ,  '  o— ' )  ; 

75  hold  on; 

76  plot  (extreme.w,  x_lower_max_global ,  '  r~  1 )  ; 

77  pause  (0.01)  ; 

78 

79  dw  =  (w.box.s i ze  (nbhd.idx) —w.box.si ze  (2  :  end)  ) /2 ; 

80 

81  if  dw(l)  >  0 

82  zero.cost  =  opt.cost.all  (2 )  ; 

83  w.box.s i ze  ( 1 )  =  w.box.size  (2 )  ; 

84  end 

85 

86  dw  =  (w.box.s i ze  (nbhd.idx) —w.box.si ze  (2  :  end)  )  ;  dw  =  dw.  *rand  (size  (dw)  )  ; 

87 

88  w.box.size  =  w_box_size  +  [0  dw  ( : )  '  ]  ; 

89  w.box.size  =  sort  (w.box.size)  ; 

90 

91  if  w.box.size  (end)  >  1 

92  w.box.size  (end)  =  1; 

93  end 

94 

95  [aa,  bb]  =min  (abs  (w.box.size— extreme.w)  )  ; 

96  w.box.si  ze  (bb)  =  extreme.w; 

97 

98  % - 

99  if  "isempty  ( w.prev.box) 

100  if  length  (w.box.size)  ==length  (w.prev.box) 

101  if  norm  (w.box.size— w.prev.box)  ==0 

102  not  .converge.sw  =  false; 

103  end 

104  end 

105  end 

106  if  length  (w.box.size )  ==1 

107  not.converge.sw  =  false; 

108  end 

109  w.prev.box  =  w_box_size; 

no  % - 

m 

112  new.w.range  =  abs  (w.box.size  (end)— w.box.size  (1)  )  ; 

113  if  abs  (new.w.range  —  current.w.range )  <  w.box.min.size 

114  count.converge  =  count.converge  +  1; 

115  else 

116  count.converge  =  0; 

117  end 

118 

119  if  new.w.range  <  w_box_min_size  ||  current.iter  >  max.iter 

120  not.converge.sw  =  false; 
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end 
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if  count.converge  >  min_iter 

w_box_size  =  extreme.w  +  extreme_w*randn  (size  (w_box_size)  )  *new_w_range/10; 
w_box_size  =  sort  ( w_box_size )  ; 

end 


%  intermediate  output 


opt_x  =  x.centre  +  extreme_w*extreme_delta  ( : )  .  *w_x  ( : )  ; 
f print f (' %d—th  iteration:  Current  Optimal  Sol:  x  =  ( ' , current-iter) ; 
for  idx=l : state.dim 

fprintf ( ' % 6 . 4 f  ' , opt_x (idx) )  ; 

end 

fprintf 

fprintf  (  ’  dw  =  %6.4f  <  %6.4f  ?\n',  abs  (w_box_size  (end)—  w_box_size  (1)  )  ,  w_box_min_si  ze ) 


current.iter  =  current.iter  +  1; 


end 
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A.  3  cal_state_bounds.m 


2  %  States  bound  calculation  for  discrete  nonlinear  systems  using  mu— bounds 

3  % - 

4  function  [x_upper,  x_lower,  extreme-delta]  =  cal_state_bounds  (k_step ,  box.init, 

5  num.steps,  num.delta.sample ,  delta.dim,  state.dim,  which.x,  ... 

6  p.centre,  w_p,  x.centre,  w_x,  safety-factor,  ... 

7  f.real.prev,  delta.px.prev,  num.replace,  ... 

8  reset.steps,  c_max,  find_min) 

9  %  Input : 

10  %  k.step:  bounds  for  k.step  value  of  f(x,d) 

n  % 

12  %  box.init :  initial  search  range  of  k,  k  is  searched  in  the  range  of 

13  %  box.init  (1)  <=  k  <=  box_init(2) 

14  % 

15  %  num.delt a.sample :  number  of  random  samples  along  the  1— norm  box  edge 

16  % 

17  %  delta.dim:  how  many  uncertain  parameters 

18  % 

19  %  state.dim:  dimension  of  state  vector  x 

20  % 

21  %  which.x:  which  x's  bounds  to  be  calculated,  range  =  [1  state.dim] 

22  % 

23  %  p.centre:  vector  for  the  centres  of  uncertain  parameters 

24  % 

25  %  w_p :  weighting  vector  for  the  uncertain  parameters  so  that  p.delta  is 

26  %  between  —1  and  1 

27  % 

28  %  x.centre:  vector  for  the  centres  of  state  vector 

29  % 

30  %  w.x :  weighting  vector  for  the  state  vector  so  that  x.delta  is 

31  %  between  —1  and  1 

32  % 

33  %  safety.f actor :  safety  factor  for  the  upper  bound,  greater  than  1, 

34  %  smaller  value  may  give  more  tighter  bound  but  with  higher  risk 

35  %  that  may  be  wrong,  this  only  applies  to  the  upper  bound 

36  % 

37  %  c.min.max :  adjusting  value  for  obtain  bounds  for  max(f)  and  min(f) 

38  % 

39  %  Output : 

40  %  x.upper:  upper  bound 

41  %  x.lower:  lower  bound 

42  %  f.min.max:  minimum  and  maximum  values  of  f  found  during  the  sampling 

43 

44  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

45  %  University  of  Glasgow 

46  %  Biomedical  Engineering 

47  %  RobustLab:  http://www.robustlab.org 

48 

49  %% 

50  tol  =  le— 3; 

51 

52  %% 

53  current.big.box  =  box.init  (2)  ; 

54  current.small.box  =  box.init  (1)  ; 

55  f_min_max  =  [inf  —inf]  ; 

56  current.delt  a.px.wor  st  =  inf; 

57  x.lower  =  —inf; 

58 
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if  isempty (delta.dim) 
delta-dim  =  0; 

end 

f  easible_set_sw  =  true; 

while  isinf  ( current_delta_px_worst  (1)  ) 

init-box_saved  =  [ current_small_box  current_big_box]  ; 

while  abs  ( current_big_box— current_small_box)  >  tol 

current_box_size  =  ( current_big_box  +  current_small_box)  /2 ; 
current_k  =  current_box_size; 

%  uniform  sample  in  the  1— norm  box  random  face 
delta.x  =  2*rand  (num_delta_sample ,  delta.dim+state.dim,  1)  —  1; 
rand_face  =  randi  (delta.dim+state.dim,  num_delta_sample ,  1 )  ; 
idx-face  =  sub2ind  (  [num_delta_sample  delta_dim+state_dim]  ,  ... 

( 1 :  num.de  1 1 a_s ampl e )  ' ,  rand.f  ace )  ; 
delta.x  ( idx.f  ace )  =  sign  (rand  (num_delta_sample,  1)  —  0 . 5)  ; 

k_delta_c  =  current.k; 

[f_real_org,  f.delta]  =  I.ND  (delta.x,  k.delta.c,  x. centre,  w.x,  c_max,  find_min) 
[f.real,  idx.f  _real_min  ]  =  min  (  f  _real_org)  ; 

if  f.real  <  0 

current  _delta.px.wor  st  =  delta.x  ( idx.f  .real.min ,  :)  ; 

current  _big_box  =  current  Jsox.size; 

cand.k  =  1/ (abs  ( f.delta  ( idx.f  .real.min  ,  :)))  ; 
if  cand.k  <  current.k 
current.k  =  cand.k; 
x.lower.cand  =  l/cand_k; 

else 

x.lower.cand  =  l/current_k; 
current.k.upper  =  current.k; 

end 

if  x.lower.cand  >  x.lower 
x.lower  =  x.lower.cand; 
current.k.upper  =  current.k; 
extreme.delta  =  delta.x  ( idx.f _real_min ,  :)  ; 

end 

end 

if  f.real  >  0 

current.small.box  =  current  Jsox.size; 

end 

%  no  feasible  set 

if  sum (abs ( f.delta) ) ==0  &&  length ( f.delta)  >  1 
x.upper  =  0; 
x.lower  =  0; 

extreme.delta  =  delta.x (1, :) ;  %  arbitrary  delta  for  flat  cost  function 

f easible.set.sw  =  false; 

break; 

end 

end  %  of  while:  check  if  the  boxes  are  close  enough? 
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%  if  never  find  f_real  <  0,  then  adjust  the  initial  box  sizes 
current_big_box  =  init_box_saved  (2 )  *2; 
current_small_box  =  max  ([1/pi  init_box_saved  ( 1 ) /2  ]  )  ; 

if  ~ f easible_set_sw 
break; 

end 

end  %  of  current_delta_px_worst  ==  inf? 

%  Refining  the  x.upper,  i.e.,  k.lower 
if  f easible_set_sw 

upper_bd_not  .converge  =  true; 

epsilon_box  =  tol*10; 

k_box  =  [ current_k_upper/10  current_k_upper]  ; 
while  upper_bd_not -converge 
for  idx=l :  num.steps 

%  sample  in  the  epsilon— norm  box  random  face 

delta_x  =  2*rand  (num_delta_sample*10,  delta_dim+state_dim,  1)  —  1; 
delta.x  =  epsilon_box*delta_x; 

delta_x  =  kron  (ones  (num_delta_sample*10, 1)  ,  current_delta_px_worst ) 
k_delta_c  =  k_box(l)  +  abs  (dif f  (k_box)  ) /2 ; 

[f_real.org,  f.delta]  =  I_ND  (delta.x,  k_delta_c,  x.centre,  ... 
w_x,  cjnax,  find_min)  ; 

f_real  =  min  (  f  _real_org)  ; 
f_delta_min  =  min  ( f.delta,  [ ]  ,  1 )  ; 
f_delta_max  =  max  ( f_delta,  [ ]  ,  1 )  ; 

if  f_delta_min  <  f _min_max  ( 1 ) 

f_min_max(l)  =  f_delta_min; 

end 

if  f_delta_max  >  f_min_max(2) 

f_min_max(2)  =  f_delta_max; 

end 

if  f.real  <  0 

k_box(2)  =  k_delta_c; 
break; 

end 


end 

if  (  (idx==num_steps)  &&  (f.real  >  0)  ) 
upper_bd_not_converge  =  false; 
x.upper  =  l/k_box(2); 

end 

end  %  of  while 

end  %  of  if  ~ f easible.set.sw 
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delta_x 


A. 4  IJMD.m 


1  function  [f.real,  f_delta]  =  I.ND  (delta.x,  k.delta.c,  x.centre,  w.x,  c.max,  find_min) 

2  %  Input : 

3  %  k.step:  bounds  for  k.step  value  of  f(x,d) 

4  % 

5  %  delta_px:  vector  including  delta_p  and  delta_x 

6  %  e.g.  delta.px  =  [  delta_p(l);  delta_p(2);  delta.xl;  delta_x2] 

7  % 

8  %  k_delta_c:  current  k*delta_c 

9  % 

10  %  which.x:  which  x's  bounds  to  be  calculated,  range  =  [1  state_dim] 

n  % 

12  %  p.centre:  vector  for  the  centres  of  uncertain  parameters 

13  % 

14  %  w_p :  weighting  vector  for  the  uncertain  parameters  so  that  p_delta  is 

15  %  between  —1  and  1 

16  % 

17  %  x.centre:  vector  for  the  centres  of  state  vector 

18  % 

19  %  w.x :  weighting  vector  for  the  state  vector  so  that  x.delta  is 

20  %  between  —1  and  1 

21  % 

22  %  c.min.max:  adjusting  value  for  obtain  bounds  for  max(f)  and  min(f) 

23  % 

24  %  Output : 

25  %  f  .real :  real  part  of  det(I— N  Delta) 

26  %  f .delta:  f  evaluated  at  the  sampled  point 

27 

28  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

29  %  University  of  Glasgow 

30  %  Biomedical  Engineering 

31  %  RobustLab:  http://www.robustlab.org 

32 

33  f .delta  =  get.f .delta  (delta.x,  x.centre,  w.x,  c_max)  ; 

34  f.real  =  1  —  k_delta_c*abs  ( f .delta)  ; 

35 

36  end 
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A. 5  get_f_delta.m 


1  function  f.delta  =  get_f _delta  (delta.x,  x.centre,  w_x,  c_max) 

2  %  Input : 

3  %  k.step:  bounds  for  k.step  value  of  f(x,d) 

4  % 

5  %  delta_px:  vector  including  delta_p  and  delta_x 

6  %  e.g.  delta.px  =  [  delta_p(l);  delta_p(2);  delta.xl;  delta_x2] 

7  % 

8  %  p.centre:  vector  for  the  centres  of  uncertain  parameters 

9  % 

10  %  w_p :  weighting  vector  for  the  uncertain  parameters  so  that  p_delta  is 

n  %  between  —1  and  1 

12  % 

13  %  x.centre:  vector  for  the  centres  of  state  vector 

14  % 

15  %  w_x :  weighting  vector  for  the  state  vector  so  that  x.delta  is 

16  %  between  —1  and  1 

17  % 

is  %  Output : 

19  %  f_delta:  function  f  evaluated  at  the  sampled  x  and  p 

20 

21  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

22  %  University  of  Glasgow 

23  %  Biomedical  Engineering 

24  %  RobustLab:  http://www.robustlab.org 

25 

26  %%  DO  NOT  CHANGE  HERE 

27 

28  %  pertubed  states 

29  n.delta  =  size  (delta.x,  1 )  ; 

30  delta_dim  =  size  (delta_x,  2 )  ; 

31 

32  x.centre  =  x.centre  ('•)'? 

33  w_x  =  w_x  ( :  )  ' ; 

34  x_delta  =  kron  (ones  (n.delta,  1)  ,  x.centre)  +kron  (ones  (n.delta,  1)  ,  w.x)  .  *delta_x 

35 

36  %%  ONLY  CHANGE  BELOW 

37  % - 

38  %  f  (X,  p) 

39  % - 

40  f_delta  =  zeros  (size  (x.delta,  1)  ,  1)  ; 

41 

42  for  idx=l :  1 :  delta_dim— 1 

43  xl  =  x_delta  ( : ,  idx)  ; 

44  x2  =  x.delt a  ( : ,  idx+1 )  ; 

45  f -delta  =  f -delta  +  100*  (x2-xl .  ~2 )  .  "2  +  (l-xl).~2; 

46  end 

47 

48  f_delta  =  abs(f_delta  +  c_max)  ; 

49  % - 

50 

51  end 
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Appendix  B 

ErbB  Signalling  Pathways: 
GPU  Example 

B.l  state  bounds  for  many  steps  main. m 


2  %  Main:  States  bound  calculation  for  many  k  steps 

3  % - 

4  % 

5  %  See  #1,  #2  and  #3  below 

6 

7  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

8  %  University  of  Glasgow 

9  %  Biomedical  Engineering 

10  %  RobustLab:  http://www.robustlab.org 

n 

12  clear  all; 

13 

14  global  f unct ion_handle ; 

15 

16  num_data  =  1024;%*10; 

17  funct ion_handle=  .  .  . 

is  parallel .  gpu  .  CUDAKernel  (  '  get.f_delta_ptx_vector.form  .  ptx ' ,  ’  get.f_delta_ptx_vector.form  .  cu  ' 

19  f unction_handle  .  GridSize  ( 1 )  =  ceil  (num_data/function_handle  . MaxThreadsPerBlock)  ; 

20  num.thread  =  min  (  [num.data  function.handle  .MaxThreadsPerBlock]  )  ; 

21  funct ion_handle  .  ThreadBlockSize  ( 1 )  =  num.thread; 

22 

23  %% 

24  % - 

25  %  1 .  Algorithm  Parameters  that  you  may  want  to  adjust 

26  % - 

27  load  init ial.condit ion; 

28 

29  dt=0.1; 

30 

31  %  number  of  samples  in  one  Period 

32  num.steps  =  1; 

33  num.delta.sample  =  num.data;  %1024*10; 

34  %  lower  with  less  computation  but  higher  risk  to  be  failed 

35  %  higher  with  less  risk  but  more  computational  burden 

36 
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38 


saf  ety_f  actor_min 


39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
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82 
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84 

85 
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87 

88 

89 
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91 

92 

93 

94 

95 

96 

97 

98 


=  1.00  +  le— 9; 

always  greater  than  or  equal  to  1 
smaller  would  give  tighter  upper  bounds 
failed 

larger  would  have  less  possbility  to  be 
tight  upper  bounds 


with  higher  risk  to  be 
failed  but  less 


s.cale  =  2; 

%  uncertain  parameters  dimension 

delta_dim  =  227;  %  how  many  uncertain  parameters 
state_dim  =  504;  %  dimension  of  state  vector  x 
which_x  =1;  %  which  x’s  bounds  to  be  calculated 

%  Range  of  k,  where  k  for  k*delta_c 
box.init  =  [le— 20  le20] ; 

%  you  should  know  roughy  where  the  k.optimal  would  be. 
%  i.e,  here  I  am  quite  sure  that  le— 2  <  f  <  lel2 
%  if  your  calculation  failed,  then  you  may  adjust  this 


%  2.  the  function  f(x,p)  related  setting 


%  what  k— th  step  value  bounds  for  x_k  =  f(x_k— 1,  p) 
k_step_all  =  1:300; 
reset_steps  =  3000; 

%  parameter  bounds : 

get_p_nominal ;  %  p_centre  from  here 
w_p  =  p_centre*0 . 1 ; 

%  state  bounds 

w_x  =  (  (double  (x_initial==0 )  +x_initial)  *le—  6)  *0; 
w_x  ( w_x>le— 6)  =le— 6*0; 
x.centre  =  x.initial; 

w.x.prev  =  w_x; 

x_cent  re.prev  =  x.centre; 

w_x.org  =  w_x; 
x_centre.org  =  x.centre; 


%  3.  Finally,  you  may  want  to  use  your  own  f(x,p)  function. 
%  Then,  update  "get.f .delta . m" 


%  No  need  to  change  from  here 


for  k.step  =  k.step.all 


%  reset  uncertain  parameter  range 


if  (mod  (k.step,  reset.steps)  ==1)  &&  (k.step  >  1) 

x.centre.t  (which.x)  =  (x_upper_max_all  (end)  +  x.upper.min.all  (end)  ) /2 ; 
w.x.t  (which.x)  =  ( x_upper.max.al  1  (end)  —  x.upper.min.all  (end)  ) /2; 
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99 


x_remain=setdif  f  (1 :  state.dim,  which.x)  ; 
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for  x_step=l : length (x.remain) 

current.x  =  x_remain  (x_step)  ; 

[x.upper.t ,  tl,  t2]  =  cal_state_bounds_gpu  (reset.steps,  ... 

box.init ,  num.steps ,  num.delta.s ample ,  delta.dim,  state_dim,  .  .  . 
current.x,  p.centre ,  w_p ,  x.cent re.prev,  w.x.prev ,  .  .  . 

safety.f actor , reset.steps *2 , [ ] , dt) ; 

c_max_t  =  2  *  x.uppe  r.t ;  c.min.t  =  — 2*x_upper_t ; 

[x_upper_max_t ,  tl,  t2]  =  cal_state_bounds_gpu  (reset.steps ,  ... 

box_init ,  num.steps ,  num_delt a.sample ,  delt a.dim,  state_dim,  .  .  . 
current.x,  p.centre ,  w_p ,  x.cent  re.prev,  w.x.prev ,  .  .  . 

safety-factor,  reset_steps*2 ,  c_max_t ,  dt)  ; 

x.uppe  r  .max.t  =  x.uppe  r  .max.t  —  c.max.t; 

[x.upper.min.t ,  tl,  t2  ]  =  cal.state_bounds.gpu  (reset.steps ,  ... 

box.init ,  num.steps ,  num.delt a.sample ,  delt a.dim,  state.dim,  .  .  . 
current.x,  p.centre ,  w_p ,  x.cent  re.prev,  w.x.prev ,  .  .  . 

safety.f actor ,  reset_steps*2 ,  c.min.t ,  dt)  ; 

x_upper_min_t  =  —  x_upper_min_t  —  c.min.t ; 

x.cent re_t  ( current.x)  =  (x.upper.max.t +  x_upper_min_t )  /2; 
w.x.t  (current.x)  =  (x.upper.max.t— x.upper.min.t )  /2 ; 


end 

x.centre  =  x.cent re_t  (:)  ; 
x.cent  re.prev  =  x.centre; 
w_x  =  w.x.t  ( : )  ; 
w.x.prev  =  w.x; 


end 


%  Calculate  Bounds 


%  calculate  c_max  and  c_min 


fprintf  (  ' - - -  -  \  i i '  )  ; 

fprintf ( ’  Calculate  c.max  and  c_min\n'); 

fprintf ( ' ==============================================\n ' )  ; 


[x.upper,  x.lower,  f_min_max_l]  =  cal.state.bounds.gpu  (k.step ,  ... 

box.init ,  num.steps ,  num.delta.sample ,  delta.dim,  state.dim,  .  .  . 
which.x ,  p.centre ,  w_p ,  x.centre ,  w.x ,  safety-factor,  reset.steps ,  [  ]  ,  dt )  ; 
c.max  =  s_scale*x_upper;  %  to  just  make  sure  all  positive  so  it's  doubled 


c.min  =  -s_scale*x_upper;  %  for  the  similar  reason,  make  it  all  negative 

%  calcualte  upper  and  lower  bounds  for  the  maximum 

fprintf ( ' ==============================================\n ' )  ; 

fprintf ( ’  Calculate  bounds  for  maximum\n ' ) ; 

fprintf ( ' ==============================================\n ’ )  ; 

[  x_upper_max ,  x_lower_max,  f_min_max_2]  =  cal.state.bounds.gpu  (k.step,  ... 
box.init ,  num.steps ,  num.delta.sample ,  delta.dim,  state.dim,  .  .  . 
which.x ,  p.centre ,  w.p ,  x.centre ,  w.x ,  saf ety.f actor  jnax ,  reset.steps ,  c_max ,  dt ) 
x.uppe  r  .max  =  x.uppe  r  .max  —  c.max ; 
x_lower_max  =  x_lower_max  —  c.max; 

%  calculate  upper  and  lower  bounds  for  the  minimum 

fprintf ( ' ==============================================\n  '  )  ; 

fprintf ( 1  Calculate  bounds  for  minimum\n ' ) ; 
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fprintf ( ' ==============================================\n ' )  ; 

[x_upper_min,  x_lower_min,  f_min_max_3]  =  cal.state_bounds.gpu  (k.step,  ... 
box.init ,  num.steps ,  num.de It a.s ample ,  delta.dim,  state  .dim,  .  .  . 
which.x ,  p.centre ,  w_p ,  x.centre ,  w_x ,  saf ety_f actorjnin ,  reset.steps ,  cjnin ,  dt )  ; 
x_upper_min  =  — x_upper_min  —  c_min; 
x.lower.min  =  —  x.lower  jnin  —  c.min; 

f  .min.max.by.simulat ion  =  [  f  .min.max.l ;  f_min_max_2;  f  _min_max_3  ]  ; 
f  .min.max.by.simulat  ion  =  ... 

[min  ( f_min_max_by_s imulat  ion  ( : ,  1 )  )  max  ( f_min_max_by_s imulat ion  ( : ,  2 )  )  ]  ; 

x_upper_max_all  (k_step)  =  x_upper_max; 
x_lower.max.all  (k.step)  =  x_lower_max; 

x_upper_min_all  (k.step)  =  x_upper_min; 
x_lower_min_all  (k_step)  =  x_lower_min; 

end 

figure ; 

%  maximum 

k_step_all  =  [0  k_step_all]; 

plot  (k_step_all*dt,  [x .centre _org  (which.x)  +w_x_org  (which.x)  x_upper_max_all  ]  ,  '  ro— '  ) 
hold  on; 

%  minimum 

plot  (k_step_all*dt,  [x_centre.org  (which.x)-w.x.org  (which.x)  x_upper_min_all]  ,  ' bx— ') 
xlabel ( ' k ' ) ; 

ylabel_str  =  sprintf (' State  x%d' ,  which.x)  ; 
ylabel (ylabel.str)  ; 
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B.2 


cal_state_bounds_gpu.m 


2  %  States  bound  calculation  for  discrete  nonlinear  systems  using  mu— bounds 

3  % - 

4  function  [x_upper,  x_lower,  f_min_max]  =  cal_state_bounds_gpu  (k.step ,  box.init, 

5  num.steps,  num.delta.sample ,  delta.dim,  state.dim,  which.x,  ... 

6  p.centre,  w_p  ,  x.centre,  w.x,  saf  ety_f  actor ,  re  set  .steps,  c_min_max,  dt) 

7  %  Input : 

8  %  k.step:  bounds  for  k.step  value  of  f (x, d) 

9  % 

10  %  box.init :  initial  search  range  of  k,  k  is  searched  in  the  range  of 

n  %  box.init  (1)  <=  k  <=  box.init  (2) 

12  % 

13  %  num.delta.sample:  number  of  random  samples  along  the  1— norm  box  edge 

14  % 

15  %  delta.dim:  how  many  uncertain  parameters 

16  % 

17  %  state.dim:  dimension  of  state  vector  x 

18  % 

19  %  which.x:  which  x's  bounds  to  be  calculated,  range  =  [1  state.dim] 

20  % 

21  %  p.centre:  vector  for  the  centres  of  uncertain  parameters 

22  % 

23  %  w.p :  weighting  vector  for  the  uncertain  parameters  so  that  p.delta  is 

24  %  between  —1  and  1 

25  % 

26  %  x.centre:  vector  for  the  centres  of  state  vector 

27  % 

28  %  w.x :  weighting  vector  for  the  state  vector  so  that  x.delta  is 

29  %  between  —1  and  1 

30  % 

31  %  saf ety.f actor :  safety  factor  for  the  upper  bound,  greater  than  1, 

32  %  smaller  value  may  give  more  tighter  bound  but  with  higher  risk 

33  %  that  may  be  wrong,  this  only  applies  to  the  upper  bound 

34  % 

35  %  c.min.max:  adjusting  value  for  obtain  bounds  for  max(f)  and  min(f) 

36  % 

37  %  Output : 

38  %  x.upper:  upper  bound 

39  %  x.lower:  lower  bound 

40  %  f .min.max :  minimum  and  maximum  values  of  f  found  during  the  sampling 

41 

42  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

43  %  University  of  Glasgow 

44  %  Biomedical  Engineering 

45  %  RobustLab:  http://www.robustlab.org 

46 

47  %% 

48  if  nargin  <  13 

49  reset.steps  =  1000; 

50  end 

51  if  nargin  <  14 

52  c.min.max  =  0; 

53  end 

54  if  nargin  <  15 

55  dt  =  1 ; 

56  end 

57 

58  if  isempty  ( c.min.max ) 
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59 


c_min_max  =  0; 


60  end 

61 

62  tol  =  le— 6; 

63  current_big_box  =  box_init  (2)  ; 

64  current _small_box  =  box.init  (1)  ; 

65  f_min_max  =  [inf  —inf]  ; 

66 

67  %% 

68  while  abs  ( current_big_box— current_small_box)  >  tol 

69 

70  current_box_size  =  ( current_big_box  +  current_small_box) /2 ; 

71  current_k  =  current_box_size; 

72 

73  comp.pct  =  abs  (current_big_box— current_small_box)  ; 

74 

75  if  comp.pct  <10  &&  comp.pct  >  le— 3 

76  fprintf ( ' (%dth— step)  Remaining  calulation  =  %4.3f%%  (a . u . ) \n ' , k.step , comp.pct ) 

77  end 

78 

79  for  idx=l :  num_steps 

80 

81  %  sample  in  the  1— norm  box  random  face 

82  delta_px  =  2*rand  (num_delta_sample ,  delta_dim+state_dim)  —  1; 

83  rand_face  =  randi  (delta_dim+state_dim,  num.delta.sample ,  1)  ; 

84  idx.face  =  ... 

85  sub2ind  (  [num.delta.sample  delta.dim+state.dim]  ,  ( 1 :  num.delta.sample )  ' ,  rand.face)  ; 

86 

87  delta.px  ( idx.f  ace )  =  sign  (rand  (num.delta.sample,  1)  —  0 . 5)  ; 

88  k_delta_c  =  current_k; 

89 

90  [f.real,  f_delta]  =  I_ND_gpu_vector  (k_step ,  delta.px,  k_delta_c ,  ... 

91  which_x ,  p.centre  ,w.p,  x.centre ,  w_x ,  c_min_max ,  re  set  .steps ,  delta.dim,  dt )  ; 

92 

93  [f.real,  idx_f  _real_min  ]  =  min(f.real); 

94  f_delta_min  =  min  ( f .delta ,  [],  1 )  ; 

95  f_delta_max  =  max  ( f .delta ,  [],  1 )  ; 

96 

97  if  f_delta_min  (which.x)  <  f _min_max  (1 ) 

98  f.min.max(l)  =  f .delta  (which.x)  ; 

99  end 
100 

101  if  f.delta.max  (which.x)  >  f_min_max(2) 

102  f_min_max(2)  =  f .delta  (which.x)  ; 

103  end 

104 

105 

106  if  f.real  <  0 

107  current  _big_box  =  current  Jsox.size; 

108  x.lower  =  1/ current.k; 

109  current  .k.upper  =  current.k; 

no  current _delta.px.wor st  =  delta.px  ( idx.f  .real.min ,  : )  ; 

m  break; 

112  end 

113 

114  end 

115 

116  if  (idx==num_steps )  &&  (f.real  >  0) 

117  current.small  Jdox  =  current_box_size; 

ns  x.upper  =  1/ current.k; 

119  end 

120 


66 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited 


121 
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%  Refining  the  x.upper,  i.e.,  k.lower 

ext reme_delta  =  sign  (current_delta_px_worst )  ; 

upper_bd_not_converge  =  true; 

epsilon_box  =  min([tol*10  norm  (extreme_delta— current.delta.px.worst )  ]  )  ; 
k_box  =  [  current  .k.upper /10  current.k.upper] ; 

while  upper_bd_not .converge 


for  idx=l : num_steps 

%  sample  in  the  epsilon— norm  box  random  face 

delta.px  =  2*rand  (num_delta_sample ,  delta.dim+state.dim,  1)  —  1; 
delta.px  =  epsilon_box*delta_px; 

delta.px  =  kron  (ones  (num.delta.sample,  1)  ,  current.delta.px.worst)  +  delta.px; 

k_delta_c  =  k_box(l)  +  abs (dif f (k_box) ) /2 ; 

[f.real,  f_delta]  =  I_ND_gpu_vector  (k_step ,  delta.px,  ... 

k_delt  a_c ,  which.x ,  p.centre ,  w_p ,  x.centre ,  w.x ,  c_min_max ,  re  set  .steps ,  delt  a.dim,  dt ) 

f.real  =  min(f.real); 
f_delta_min  =  min  ( f .delta ,  [],  1 )  ; 
f_delta_max  =  max  ( f .delta ,[],  1 )  ; 

if  f.delta.min  (which.x)  <  f_min_max(l) 
f.min.max(l)  =  f.delta  (which.x)  ; 

end 

if  f.delta.max  (which.x)  >  f_min_max(2) 
f_min_max(2)  =  f.delta  (which.x)  ; 

end 

if  f.real  <  0 

k_box(2)  =  k.delta.c; 
break; 

end 

end 

if  ( (idx==num_steps)  &&  (f.real  >0))%  ||  (abs (dif f (k.box) )  <  tol) 
upper_bd_not  .converge  =  false; 
x.upper  =  l/k_box(2); 

end 


end  %  of  while 


x.upper  =  safety.f  actor  *x_upper ; 
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B.3  I_ND_gpu_vector.m 


1  function  [f.real,  f_delta]  =  I_ND_gpu_vector  (k.step ,  delta_px,  k_delta_c,  ... 

2  which_x,  p_centre,  w_p,  x_centre,  w_x,  c_min_max,  reset.steps,  num_delta ,  dt ) 

3  %  Input : 

4  %  k_step:  bounds  for  k_step  value  of  f(x,d) 

5  % 

6  %  delta.px:  vector  including  delta.p  and  delta.x 

7  %  e.g.  delta.px  =  [  delta_p(l);  delta_p(2);  delta.xl;  delta_x2] 

8  % 

9  %  k_delta_c:  current  k*delta_c 

10  % 

n  %  which_x:  which  x's  bounds  to  be  calculated,  range  =  [1  state_dim] 

12  % 

13  %  p.centre:  vector  for  the  centres  of  uncertain  parameters 

14  % 

15  %  w_p :  weighting  vector  for  the  uncertain  parameters  so  that  p.delta  is 

16  %  between  —1  and  1 

17  % 

is  %  x_centre:  vector  for  the  centres  of  state  vector 

19  % 

20  %  w_x :  weighting  vector  for  the  state  vector  so  that  x_delta  is 

21  %  between  —1  and  1 

22  % 

23  %  c_min_max :  adjusting  value  for  obtain  bounds  for  max(f)  and  min(f) 

24  % 

25  %  Output : 

26  %  f .real :  real  part  of  det(I— N  Delta) 

27  %  f_delta:  f  evaluated  at  the  sampled  point 

28 

29  %  Programmed  by  Jongrae  <Jongrae . Kim@glasgow . ac . uk> 

30  %  University  of  Glasgow 

31  %  Biomedical  Engineering 

32  %  RobustLab:  http://www.robustlab.org 

33 

34  global  f unct ion_handle ; 

35 

36  k_step  =  mod  (k_step,  reset_steps)  ; 

37  if  k_st ep==0 

38  k_step=reset_steps ; 

39  end 

40 

41  %%  CHANGE  BELOW 

42  % - 

43  num.state  =  length  (x_centre )  ; 

44  num.data  =  size  (delta_px,  1 )  ; 

45 

46  %  output  state 

47  switch  which.x 

48  case  1 


49 
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50 

92 

93 

98 

99  104  105  110 
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m 

116 

117 

122 

123 

128 

129 

134 

135 

212 

215 

218  221 

224 

227 

230 

233 

264 

265 

280 

283 

296 

299 

309 

312 

377 

378 

390 

410 

436 

446 

447 

459 

465 

477  ] 

52 

ErbB14  =  [ 

51 

52 

94 

95 

100 

101 

106 

107 

112 

113 

118 

119 

124 

125 

130 

131 

136 

137 

211 

214 

217  220 

223 

226 

229 

232 

266 

267 

281 

284 

297 

300 

310 

313 

373 

374 

391 

411 

437 

448 

449 

458 

466 

478] 

53 

out. 

_idx  = 

[ErbBl 

ErbBl 

ErbB12  ErbB13  : 

ErbBl 4 

] ;  %  ErbBl  is 

dimer 

54 

case  2 

55 

Erk 

=  [346 

i  346 

348 

349 

359 

360 

361 

362  363  364 

442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453  454 

455]  ; 

56  out.idx  =  Erk; 

57  case  3 

58  Akt  =  [405  407  470  471]; 

59  out_idx  =  Akt; 

60  otherwise 

61  out.idx  =  []; 

62  end 

63 

64  %  if  the  centre  of  x  is  zero,  then  take  the  absolute  value  of  all 

65  %  perturbation 

66 

67  big_w_p  =  kron  (ones  (num.data,  1)  ,  w_p  ( : )  '  )  ; 

68  %big_p  =  kron  (ones  (num.data,  1)  ,  p.centre  (:)' )  + 

69  big_p  =  big_w_p  .  *delt a_px  ( : ,  1 :  num.delta )  ; 

70  big_w_x  =  kron  (ones  (num.data,  1)  ,  w.x  ( : )  '  )  ; 

71  big.x.c  =  kron  (ones  (num.data,  1)  ,  x_centre  (:)'); 

72 

73  %  take  absolute  perturbation  if  x_c  =  0 

74  zero.idx  =  (big_x_c==0)  *  (  —  1)  +  (big_x_c~=0)  ; 

75  x.state  =  big.x.c  +  big.w.x  .  *  (delta.px  ( :  ,  num_delta+l :  end)  .  *zero_idx)  ; 

76  x.state  (x_state<0)  =  0; 

77 

78  %keyboard 

79 

so  %  include  nominal  for  safety 
si  x.state  (end,  : )  =  x.cent re  ( : )  ' ; 

82  big.p  (end,  : )  =big_p  (end,  : )  *0; 

83 

84  p.gpu  =  gpuArray  (big.p  (:))  ; 

85  x.state.gpu  =  gpuArray  (x.state  (:))  ; 

86  k.step.gpu  =  gpuArray  (k.step)  ; 

87  num.data.gpu  =  gpuArray  (num.data)  ; 

88 

89  %dt  =  dt* 0.001; 

90  dt.gpu  =  gpuArray (dt ) ; 

91 

92  x.out  =  zeros  (size  (x.state ))  ; 

93  x.gpu.out  =  gpuArray  (x.out )  ; 

94 

95  fO  =  zeros  (size  (x.state ))  ; 

96  fO.gpu  =  gpuArray  (fO )  ; 

97  fl  =  zeros  (size  (x.state ))  ; 

98  fl.gpu  =  gpuArray (fl) ; ; 

99  xO  =  zeros  (size  (x.state ))  ; 

100  xO.gpu  =  gpuArray  (xO )  ; 

101  xl  =  zeros  (size  (x.state ))  ; 

102  xl.gpu  =  gpuArray  (xl )  ; 

103 

104  % - 

105  %  ptx  option:  !nvcc  — ptx  get.f_delta_ptx.cu 
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106 


107  x_gpu_out  =  ... 

108  f  eval  ( f  unct  ion.handle ,  ... 

109  x.gpu.out,  p.gpu,  x.state.gpu,  k.step.gpu,  ... 

no  num_data_gpu,  dt_gpu,  fO.gpu,  fl.gpu,  xO.gpu,  xl_gpu)  ; 

m  x.state  =  gather  (x.gpu.out )  ; 

112  x_state  =  reshape  (x.state ,  num.data,  num_state)  ; 

113  f  .delta  =  sum  (x.state  (:,  out.idx)  ,  2 )  ; 

114 

115  x_state_gpu  =  gpuArray  (x.state  (:))  ; 

116 

117  f.real  =  1— k_delta_c*abs  ( f .delta  (:)  +c_min_max )  ; 

118 

119 

120  end 
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B.4  get_f_delta_ptx_vector_form.cu 


1  //  Programmed  by  Jongrae  Kim  <Jongrae . Kim@glasgow . ac . uk> 

2 

3  //  Warning 

4  //  1.  Do  not  use  a  variable  as  an  input  and  output  at  the  same  time 


#def ine  NUM.STATE 
#def ine  NUM.DELTA 

//  function  header 


504 

227 


10 

_ device _ void  f.p.x  (double 

*  x.out , 

11 

const  double  *  p_ 

in,  const 

double 

12 

13 

//  main  function 

14 

..global.,  void  get.f.delta  ( 

double 

*  x2 , 

15 

const 

double  * 

p.in,  ' 

16  const  int  k.step,  const  int  num_data,  const  double  dt, 

17  double  *  fO,  double  *  fl,  double  *  xO,  double  *  xl) 

IS  { 

19  int  idx  =  threadldx.x  +  blockldx.x  *  blockDim.x; 

20  int  currents ime_step ; 

21  int  current.state.step; 


26 

27 


//  integration  using  Two-Step  Adams— Bashforth 

for  ( current.t ime_step=0 ;  current.t ime.step  <  k.step;  current.t ime_step++) 

{ 

if  ( current_time_step  ==  0) 

{ 

//  fO  =  f(xO) 

f.p.x  (fO,  p.in,  x_in,  num.data,  idx); 


33 

34 


//  x2  =  xO  +  dt  f (xO) :  special  case  for  current_time_step  =  0 

for  (current_state_step=0;  current.state.step  <  NUM.STATE;  current_state_step++)  { 
x2  [idx+current_state_step*num_data]  =  x_in  [  idx+current_state_step*num_data  ] 

+  dt*f0  [idx+current_state_step*num_data]  ; 


37 

38 


if  (x2  [idx+current_state_step*num_data]  <  0.0) 
x2  [idx+current_state.step*num.data]  =  0.0; 


47 

48 


} 

} 

else 

{ 

ii- 


ii- 


//save  for  the  next  step 

xO  [idx+current_state_step*num_data]  =  x_in  [ idx+current_state_step*num_data ]  ; 
xl  [idx+current_state_step*num_data]  =  x2  [idx+current_state_step*num_data]  ; 


//  Two— step  Adams— Bashforth  Integration 


51 

52 


//  f(x0)  &  f ( x 1 ) 

f.p.x  (fO,  p.in,  xO,  num.data,  idx); 
f.p.x (fl,  p.in,  xl,  num.data,  idx); 


//  x_2  =  xl  +  1.5*dt*f(xl)  -  0 . 5*dt*f (xO ) 

for  (current.state.step=0;  current.state.step  <  NUM.STATE;  current_state_step++)  { 
x2  [ idx+current_state_step*num_data ]  =  xl  [idx+current_state_step*num_data] 

+  1 . 5*dt*f  1  [idx+current_state_step*num_data] 
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59 


—  0 . 5*dt*f 0 [idx+current_state_step*num_data] ; 


60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 


if  (x2  [idx+current_state_step*num_data]  <  0.0) 
x2  [idx+current_state_step*num_data ]  =  0.0; 

//  save  for  the  next  step 

xO  [idx+current_state_step*num_data]  =  xl  [idx+current_state_step*num_data]  ; 
xl  [  idx+current_state_step*num_data  ]  =  x2  [idx+current_state_step*num_data]  ; 

} 

} 

} 

} 


//  function  definition 
..device.,  void  f_p_x  (double  *  x.out, 
const  double  *  p_in,  const 


{ 


double 


x_in,  const  int  num.data,  const  int  idx) 


//  make  the  state  index  vector 
int  current.state.step; 
int  state.idx  [NUM.STATE]  ; 

for  ( current_state_step=0 ;  current.state.step  <  NUM.STATE;  current_state_step++) 
state.idx  [current.state.step]  =  idx  +  current_state_step*num_data; 

/ /  make  the  delta  index  vector 
int  current.delta.step; 
int  delta.idx  [NUM.DELTA]  ; 

for  ( current_delta.step=0 ;  current.delta.step  <  NUM.DELTA;  current_delta_step++) 
delta.idx  [current.delta.step]  =  idx  +  current_delta_step*num_data; 


//  kinetic  parameters 

double  kdl=3 . 300000e— 02  +  p.in [delta.idx [0] ] ; 

double  klc=8 . 000000e+02  +  p.in [delta.idx [ 1 ]] ; 
double  kdlc=l . 000000e+00  +  p.in [delta.idx [2 ]] ; 
double  kdld=l . OOOOOOe— 01  +  p.in [delta.idx [3] ] ; 

/ /  omitted  . . . 

//  differential  equations 
x.out  [ state.idx  [ 0 ]  ]  =  0.0; 

//  keep  EGF  constant,  e.g.  5e— 9  =  5  nM 

//  EGF  (See  Figure  2  in  the  Reference)  or  O.Ole— 9  (Figure  2) 

x.out  [ state.idx  [  1  ]  ]  =  —  (kl*x .in  [ state.idx  [ 0 ]  ]  *x_in  [  state.idx  [1]  ]  —  kdl*x_in  [ state.idx  [2  ]  ]  ) 
x.out  [  state.idx  [  1  ]  ]  =  x.out  [  state.idx  [  1  ]  ] 

—  (k6*x_in  [state.idx  [1]  ]—  kd6*x_in  [ state.idx  [242 ]  ]  )  ; 
x.out  [  state.idx  [  1  ]  ]  =  x.out  [  state.idx  [  1  ]  ] 

—  (kl20b*x_in  [  state.idx  [482]  ]  *x_in  [  state.idx  [  1  ]  ]  — kdl2  0  *x_in  [  state.idx  [  60  ]  ]  )  ; 
x.out  [  state.idx  [  1  ]  ]  =  x.out  [  state.idx  [  1  ]  ] 

—  (kl20b*x_in  [  state.idx  [483]  ]  *x_in  [  state.idx  [  1  ]  ]  —  kdl20 *x_in  [  state.idx  [  62  ]  ]  )  ; 
x.out  [  state.idx  [  1  ]  ]  =  x.out  [  state.idx  [  1  ]  ] 

+  (kl22*x_in  [  state.idx  [417]  ]  *x_in  [  state.idx  [35]  ]  —  kdl22 *x_in  [  state.idx  [  1  ]  ]  )  ; 

x.out  [  state.idx  [2  ]  ]  =  (kl*x_in  [  state.idx  [ 0  ]  ]  *x_in  [  state.idx  [1]  ]—  kdl*x_in  [  state.idx  [2  ]  ]  )  ; 
x.out  [ state.idx  [2 ]  ]  =  x.out  [ state.idx  [2 ]  ] 

—  (k2*x_in  [ state.idx  [2 ]  ]  *x_in  [ state.idx  [8]  ]  —  kd2*x_in  [ state.idx  [  9 ]  ]  )  ; 
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121  x_out  [ state.idx  [2 ]  ]  =  x.out  [ state_idx  [2 ]  ] 

122  —  (k2*x_in  [ state.idx  [2  ]  ]  *x_in  [ state.idx  [2]  ]  —  kd2*x_in  [ state.idx  [11]  ]  )  ; 

123  x.out  [ state_idx  [2 ]  ]  =  x.out  [ state.idx  [2 ]  ] 

124  —  (k2b*x_in  [  state.idx  [ 2  ]  ]  *x_in  [  state.idx  [  14 ]  ]  — kd2b*x_in  [  state.idx  [  15 ]  ]  )  ; 

125  x.out  [ state.idx  [2 ]  ]  =  x.out  [ state.idx  [2 ]  ] 

126  —  (k2b*x_in  [ state.idx  [2 ]  ]  *x_in  [  state.idx  [16]  ]  —  kd2b*x_in  [ state.idx  [  17 ]  ]  )  ; 

127 

128  / /  omitted  .  .  . 

129 

130  x.out  [state.idx  [503]  ]  =  (kl22*x_in  [state.idx  [498]  ]  *x_in  [state.idx  [35]  ] 

131  —  kdl22*x_in  [state.idx  [503]  ]  )  ; 

132  x.out  [state.idx  [503]  ]  =  x.out  [ state.idx  [ 503 ]  ] 

133  +  (kl23h*x_in  [  state.idx  [243]  ]  *x_in  [  state.idx  [35]  ]  —  kdl23h*x_in  [  state.idx  [ 503  ]  ]  ) 

134  } 
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